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Abstract 

This  thesis  studies  how  rocks  evolve  due  to  the  coupled  effects  of  flow  and  chemical 
reaction.  The  study  was  motivated  by  various  experimental  observations,  both  in 
igneous  and  sedimentary  rocks.  In  the  first  part  of  this  thesis,  growth  of  microscopic, 
pore-scale,  features  in  sedimentary  rocks  is  theoretically  investigated.  It  is  found, 
in  agreement  with  experiments,  that  statistical  properties  of  pore-grain  interfaces 
mirror  growth  conditions.  The  shapes  of  pore-grain  intrefaces  both  influence  and  are 
influenced  by  large-scale  transport  properties  of  the  rock.  The  second  part  of  this 
thesis  employs  analytical  methods  to  study  flow  patterns  in  melt  upwelling  beneath 
mid-ocean  ridges.  It  is  shown  that  high  permeability  chcinnels  spontaneously  form, 
allowing  for  efficient  extraction  of  melt  from  the  system.  This  result  may  aid  in 
understanding  existing  geochemical  and  geological  observations.  In  the  third  part 
of  this  thesis,  I  present  a  new  3D  computer  model  that  simulates  flow  and  reaction 
through  a  porous  matrix.  The  model  is  used  to  study  and  compare  the  different 
characteristics  of  dissolution  and  deposition,  and  to  simulate  different  settings  for 
melt  upwelling  in  the  mantle. 
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Chapter  1 


Introduction 


Motivated  by  various  observations  in  sedimentary  and  igneous  rocks,  this  thesis  stud¬ 
ies  a  few  aspects  of  the  slow  process  of  rock  formation.  Such  a  study  of  formation 
of  geological  features  is  reminiscent  of  a  murder  investigation  of  an  unknown  person: 
On  the  one  hand,  the  body  (in  this  thesis,  a  rock)  is  observed  out  of  context.  It  had 
a  long  life  previous  to  this  moment,  but  only  a  few  snapshots  in  time  are  available 
to  deduce  the  full  chain  of  events.  On  the  other  hand,  there  is  a  multitude  of  data, 
only  a  small  percent  of  which  may  be  essential  in  understanding  what  is  observed, 
but  it  is  hard  to  determine  beforehand  which  parts  to  ignore.  The  studies  in  this  the¬ 
sis  propose  that  fluid  flow  and  reaction  through  a  porous  medium  strongly  influence 
the  evolution  of  rocks  and  their  resulting  properties,  and  may  explain  some  of  the 
experimentally  observed  richness. 

Understanding  coupled  flow  and  reaction  is  important  in  a  variety  of  geological 
and  industrial  settings.  Dissolution  and  precipitation  that  occur  during  brine  flow  are 
responsible,  to  a  Icirge  degree,  for  the  formation  of  sedimentary  rocks  from  the  initial 
compacted  grciins  [63].  Both  the  flow  patterns  and  the  chemistry  within  the  earths 
mantle  are  effected  by  similar  processes  of  reactive  porous  flow  (e.g.,  [10,  19,  64,  44]), 
but  here  the  fluid  is  lava  melted  from  the  grains  at  temperatures  hundreds  of  degrees 
higher  than  in  sedimentary  rocks.  Geologists,  oil  companies,  hydrologists,  companies 
concerned  by  contaminents  and  even  coffee  percolator  manufacturers  would  like  to 
understand,  and  in  the  end  quantify,  time  dependent  changes  in  geometrical  and 
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transport  properties  of  porous  media  due  to  clogging  or  corroding  processes. 

The  study  of  the  evolution  of  porous  media  during  reactive  flow  is  also  intriguing 
on  its  own  right,  being  a  relatively  basic,  yet  not  well  understood,  physical  process. 
Because  of  the  strong  non-linearities  and  the  multiple  length  and  time  scales  involved, 
coupled  flow  and  reaction  is  most  difficult  to  tackle.  Length  scales  range  between 
micrometers,  the  scale  of  a  single  pore  where  chemiccil  reaction  rates  may  be  controlled 
by  transport  and  kinetics  at  the  pore  wall,  to  tens  of  kilometers  which  can  be  the 
sccile  for  flow  through  a  sedimentciry  basin  or  the  mantle.  Time  scales  range  just  as 
widely:  Fluid  may  flow  on  time  scales  of  hours,  but  may  react  for  years  before  any 
significant  change  in  porosity  occurs.  Flow  and  reaction  through  porous  media  is  also 
one  of  the  few  physical  systems  in  which  length  and  time  scales  continuously  change 
with  time,  mzJang  it  impossible  to  define  a  unique  set  of  non-dimensional  parameters 
to  describe  flow.  For  excimple,  eis  a  porous  media  is  dissolved  by  acid,  channels  may 
form  with  chziracteristic  length  scaled  much  larger  than  the  initial  Darcy  scale,  and 
flow  rates  orders  of  magnitude  faster  them  in  the  initial  configuration  [16]. 

Structures  formed  in  rocks  would  be  relatively  uninteresting  if  pore  fluids  were 
static  and  in  complete  equilibrium  with  surrounding  rock.  Luckily,  many  natural 
systems  such  as  sedimentary  basins  or  the  upper  mantle  are  chemically  open  systems, 
with  fluids  continuously  driven  from  here  to  there  or  back.  In  this  case,  long-range 
interactions  or  extensive  disequilibrium  are  the  rule  rather  than  the  exception.  The 
dynamical  system  cem  be  studied  via  different  approaches,  ranging  from  microscopic 
studies  [92,  36,  9]  of  reactive  “particles”  in  actual  “pores”,  to  analog  network  simula¬ 
tions  [24],  to  solutions  of  macroscopic  partial  differential  equations  [14,  57].  A  delicate 
aspect  in  all  studies  is  the  bridging  across  the  scales.  Due  to  computational  limitations 
it  is  still  difficult  to  use  microscopic  models  to  study  macroscopic  behavior  that  occurs 
on  scales  much  greater  than  a  single  pore.  It  is  similarly  impossible  to  use  macro¬ 
scopic  models  to  describe  phenomena  that  occur  on  a  pore  scale.  It  is  thus  necessary 
to  use  a  priori  constitutive  laws  in  macroscopic  equations,  which  may  be  far  removed 
from  microscopic  calculations.  For  example,  solute  concentration  within  a  single  pore 
may  be  uniform,  but  great  differences  in  mineral  compositions  may  be  found  between 
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neighboring  pores  [63].  Thus,  one  may  conclude  that  although  microscopic  diffusion 
through  pore  fluids  is  most  efficient,  the  resulting  macroscopic  diffusion  coefficient 
is  much  smaller  than  naively  expected,  and  is  in  some  way  influenced  by  interpore 
interaction. 

In  this  thesis  the  different  length  scales  axe  treated  separately.  Figure  1-1  illus¬ 
trates  the  interactions  between  the  scales  of  the  different  studies  presented  in  this 
thesis. 

In  the  first  part  (Chapter  II,  also  in  [1])  I  construct  a  physical  microscopic  model 
of  pore-grain  interface  growth  in  sedimentary  rocks,  based  on  experimental  results 
[43,  53]  that  indicate  that  pore  scale  dynamics  are  controlled  by  chemical  reaction 
rates  and  not  by  transport  rates.  The  effect  of  flow  at  this  scale  is  thus  only  in  sup¬ 
plying  unequilibrated  fluid  to  pores  and  so  allowing  the  existence  of  non-equilibrium 
features,  and  not  in  supporting  gradients  as  at  the  larger  scales.  Chapter  II  proposes 
that  sedimentary  rocks,  regardless  of  their  mineralogy,  follow  a  universal  path  of  evo¬ 
lution  in  which  their  pore-scale  statistical  characteristics  continuously  change  as  they 
approach  (but  possibly  never  achieve)  local  chemical  equilibrium.  These  changing 
statistics  are  in  turn  tied  to  changes  in  permeability,  to  form  a  closed  picture  about 
the  evolving  sedimentary  rock.  It  is  interesting  that  forming  features,  of  the  order  of 
micrometers,  both  influence  and  are  influenced  by  the  large  scale  fluid  transport. 

In  the  second  part  of  this  thesis  (Chapter  III,  also  in  [3])  I  study  macroscopic 
aspects  of  flow  of  melt,  which  is  dissolving  the  surrounding  mantle  as  it  is  upwelling 
beneath  mid-ocean  ridges.  Coupled  flow  and  reaction  in  this  case  are  shown  to  be  re¬ 
sponsible  for  spontaneous  formation  of  macroscopic  channels,  of  fast  and  slow  porous 
flow,  which  span  the  whole  region  of  upwelling.  The  spontaneous  formation  of  chan¬ 
nels  may  resolve  a  long-standing  puzzle  (e.g.,  [22])  in  the  understanding  of  melt  ex¬ 
traction  from  mid-ocean  ridges.  Chapter  III  was  proceeded  by  [36],  which  presented 
initial  investigations,  both  computational  and  experimental,  and  a  description  of  the 
relevant  geochemical  processes. 

In  the  third  part  of  my  thesis  (Chapter  IV)  I  present  a  macroscopic  computer 
model  of  flow  cind  reaction  in  2D,  and  apply  it  to  both  general  aspects  of  flow  and 
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subduction  arc 
mid-ocean  ridge  volcanism 


Figure  1-1:  (A)  On  the  smallest  scales,  the  pore  scale,  this  thesis  studies  statistical 
changes  of  pore-grain  interfaces  in  sedimentary  rocks.  Results  show  that  pore-grain 
interfaces  both  quantitatively  reflect  large  scale  changes  in  the  rock  (i.e.,  amount  of 
diagenetic  alteration),  aind  effect  the  large  sccile  transport  properties.  In  this  way  scale 
(A)  is  tied  to  scale  (B).  (B)  On  the  intermediate  scale,  this  thesis  studies  the  gen¬ 
eral  effects  of  dissolution  and  precipitation  on  evolving  porous  media.  ZD  computer 
simulations  show  that  dissolution  will  cause  formation  of  preferred  high  permeabil¬ 
ity  paths,  (as  illustrated  in  this  figure  by  a  black  high  permeability  path  formed  by 
corrosive  flow,)  while  precipitation  will  diffuse  and  homogenize  any  initiedly  preferred 
paths  of  flow  (not  illustrated  here).  (C)  On  an  even  larger  scale,  coupled  flow-reaction 
systems  in  different  settings  result  in  different  geochemical  and  geological  outcomes 
[36].  This  thesis  shows,  by  means  of  linear  analysis  cind  ZD  computer  simulations, 
that  melt  upwelling  in  conditions  believed  to  describe  midocean  ridges,  will  focus  into 
high  permeability  chcuinels,  due  to  its  corrosive  effect  on  the  matrix.  On  the  other 
hand,  simulations  mimicking  intra-plate  volcanism  will  result  in  diffuse  flow,  due  to 
the  cooling  and  crystallizing  process  induced  by  the  continents,  and  a  formation  of  an 
overpressurized  region  just  below  the  crystallizing  region.  Using  these  results,  large- 
scale  geological  and  geochemical  observations  are  understood  via  interactions  on  the 
intermediate  scale  (B). 
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reaction  and  to  some  scenaxios  of  melt  extraction  from  the  mantle.  The  model  is  new 
both  in  its  ability  to  simulate  laxge  systems  more  efficiently  than  previous  models, 
and  in  its  analyticcil  macroscopic  description  of  the  deposition  process.  In  Chapter  IV 
the  two  processes  of  deposition  and  dissolution,  and  the  organization  of  the  porous 
media  due  to  each  one  of  them,  are  studied  as  opposite  processes.  The  studies  are 
performed  in  the  limit  where  the  effect  of  flow  is  maximized.  In  this  case  it  is  shown 
that  dissolution  produces  long-range  correlations  in  porosity  and  permeability  while 
deposition  produces  negative  correlations.  The  results  of  the  experiments  performed 
for  melt  migration  agree  both  with  Chapter  III  and  with  predictions  in  [36]. 

Lastly  I  will  summcirize  the  conclusions  of  this  thesis:  Coupled  flow  and  reaction 
may  be  significant  forces  in  shaping  rocks  as  they  evolve.  This  generic  physical  process 
may  be  responsible  for  the  observed  fractal  structures  on  the  pore  scale  of  sedimentary 
rocks  [38,  85,  1],  for  formation  of  channels  in  the  mantle  [3],  or  caves  in  calcite  rocks 
[16].  It  may  influence  rates  of  lava  flow,  and  in  turn  rates  of  sea-floor  accretion  or 
eruption  of  volcaxios,  and  determine  where  dikes  will  initiate  in  the  mantle  (Chapter 
IV). 

I  find  that  coupled  flow  and  reaction  are  responsible  for  changing  the  statistical 
characteristics  of  a  porous  medium,  with  dissolution  aind  deposition  having  qualita¬ 
tively  opposite  effects.  It  may  be  possible  in  the  future  to  use  the  statistics  of  the 
“geometrical  fingerprints”  to  obtain  quantitative  constraints  on  the  processes  that 
different  rocks  have  undergone.  It  is  also  clear  that  flow  and  reaction  effect  perme¬ 
ability,  in  a  way  which  needs  to  be  further  studied.  Future  studies  require  deeper 
understanding  cind  qucintification  of  the  interaction  between  the  scales. 
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Chapter  2 


Roughness  in  rocks 


Abstract 

Recent  laboratory  measurements  have  shown  that  pore  surfaces  of  most  sedimentary 
rocks  have  a  fractal  dimension  ranging  mostly  between  2.6  and  2.8.  The  lower  and 
upper  cutoffs  for  fractal  behavior  are  10“*  ^lnd  10^  fim,  respectively.  The  fractal 
dimension  increases  with  diagenetic  alteration.  To  explain  these  measurements,  we 
construct  a  physical  model  of  mineral  deposition  and  dissolution  on  a  substrate.  We 
propose  that  when  formation  dynamics  are  reaction  controlled,  the  forming  pore-grain 
interface  can  be  described  by  a  non-linear  partial  differential  equation  for  interface 
growth.  We  construct  a  discrete  particle-deposition  model  corresponding  to  these 
dynamics.  Three-dimensional  computer  simulations  of  the  model  show  that  resulting 
pore-grain  interfaces  cire  fractal,  with  a  fractal  dimension  that  increases  from  D  w  2.63 
to  J?  ~  2.84  as  the  dissolution  rate  is  increased,  in  close  agreement  with  observations. 
Additionally,  our  model  predicts  an  increcise  in  the  amplitude  of  interface  undulations 
with  dissolution  and  fractal  dimension.  We  conclude  that  geometrical  measures  of 
pore-grciin  interfaces  cire  ein  indicator  of  the  diagenetic  history  of  sedimentary  rocks, 
and  are  related  to  large  scale  changes  in  permeability. 


2.1  Introduction 

How  can  we  better  understand  the  conditions  under  which  sedimentary  rocks  form?  In 
this  paper  we  concentrate  on  statistical  measurements  and  geochemical  observations 
to  provide  us  with  new  insight  into  formation  processes.  Specifically,  we  study  how 
dissolution,  precipitation,  weathering,  erosion,  and  other  processes  that  alter  the  pore- 
space  of  sedimentary  rocks  from  its  initial  state  (i.e.,  processes  that  cause  “diagenetic 
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alteration”),  eifFect  certain  statistical  characteristic  of  these  rocks. 

As  shown  schematiccdly  in  Figure  2-1,  pore-grain  boundaxies  observed  in  sedimen- 
teiry  rocks  are  usually  quite  “rough”  and  complex.  A  growing  body  of  measurements 
[5,  30,  94,  84,  21,  38,  39,  83]  suggest  that  in  most  sandstones  cind  shales  these  pore- 
grain  interfaces  are  fractal  for  length  scales  ranging  over  4  orders  of  magnitude,  from 
approximately  10“^  to  10^  /im.  The  measured  surface  area,  S,  of  a  fractal  interface 
has  a  power-law  dependence  on  the  lateral  extent  L  of  the  interface:  i.e.,  S{L)  oc  L^, 
where  2  <  D  <  3  is  the  interface  fractal  dimension  [91].  Thus,  the  surface  area 
of  a  fractcil  interface  increases  faster  with  L  than  if  it  were  Euclidean  (D=2),  but 
slower  than  if  it  were  a  volume  filling  object  (D=3).  Fractal  dimensions  of  pore  sur¬ 
faces  in  sedimentary  rocks  are  observed  to  range  mostly  between  2.6  and  2.8,  with  D 
increasing  with  diagenetic  alteration  [39,  38,  85]. 

In  general,  a  fractal  distribution  of  features  in  space  also  indicates  spatial  power- 
law  correlations  between  them  [89].  Specifically,  the  “density-density”  correlation 
function,  designated  here  c(r),  describes  the  correlation  between  a  scalar  property  p 
at  position  vectors  r'  and  r'  -f  r, 

c(r)  =  1/V  p{t'  -1-  r)p{r')dr',  (2.1) 

where  V  is  the  sample  volume.  When  p{r')  is  a  distribution  of  solid  (p(r')  =  1)  and 
voids  (p(r')  =  0)  in  space,  then  c(r)  is  proportional  to  the  probability  of  finding  a 
solid  object  at  position  r  -f  r',  given  that  there  is  a  solid  object  at  position  r' .  For 
a  fractal  object  with  a  fractal  dimension  D,  embedded  in  the  Euclidean  dimension 
d  =  3,  this  correlation  function  scciles  like  [89] 

c(r)  ~  r^~^.  (2.2) 

If  pore  interfaces  are  indeed  fractal,  equation  2.2  implies  a  formation  process  that 
is  responsible  for  long-range  correlations  in  space;  e.g.,  crystal  growth  at  one  point 
in  the  pore  influences  growth  at  other  points.  Hence,  in  the  study  of  the  evolution 
of  rocks  one  cannot  isolate  growth  of  single  crystals  and  hope  to  fully  characterize 
the  dynamical  formation  process.  One  must  instead  consider  the  effect  of  the  growth 
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and  dissolution  of  a  large  number  of  crystals  and  their  effect  on  one  another,  in  order 
to  understand  the  statistical  surface  measurements  and  their  implications  for  rock 
formation. 

In  this  paper  we  propose  a  physical  model  of  evolution  of  pore-grain  interfaces  to 
explain  such  long-range  correlations.  Our  goal  is  to  provide  a  link  between  formation 
dynamics  in  rocks  and  measurable  statistical  properties.  In  constructing  a  model  for 
explaining  the  existence  of  rough  interfaces  in  rocks,  we  are  guided  by  three  main 
objectives: 

1.  The  model  should  describe  non-equilibrium  growth. 

2.  The  model  should  be  independent  of  mineralogy,  as  is  indicated  by  the  range  of 
length  scales  and  diversity  of  minerals  over  which  fractal  behavior  is  observed. 
Specifically,  fractal  behavior  is  observed  in  rocks  with  many  cementing  compo¬ 
nents,  from  clays  with  crysttils  as  small  as  10“^/im  to  quartz  with  up  to  lOO^m 
crystals. 

3.  The  dynamical  model  should  be  consistent  with  known  geochemical  constraints 
for  growth.  The  results  of  the  model  should  also  be  consistent  with  available 
qualitative  and  quantitative  statistical  observations  in  sedimentary  rocks,  and 
supply  an  explanation  for  the  range  of  fractal  dimensions  measured. 

Previous  attempts  to  construct  a  model  [30,  15,  94]  partially  met  the  first  two  of 
these  objectives  but  were  unable  to  meet  the  third.  In  particular,  no  model  has  yet 
explained  and  predicted  a  range  of  observed  fractal  dimensions. 

The  construction  of  our  model  follows  from  recent  experimental  studies,  (e.g.,  [65, 
53]),  that  suggest  that  most  sedimentary  rocks  form  by  reaction-controlled  kinetics. 
Kinetics  axe  reaction-controlled  when  the  rate-limiting  step  for  interface-growth  is  the 
chemical  reaction  at  the  interface  rather  than  transport  of  mineral  to  the  interface. 
We  propose  that  when  formation  dynamics  are  reaction  controlled,  the  forming  pore- 
grain  interface  can  be  described  by  the  interface  growth  equation  derived  by  Kardar, 
Parisi  and  Zhang  (KPZ)  [1986].  This  equation  describes  the  evolution  of  an  interface 
that  grows  everywhere  in  the  direction  normal  to  the  interface,  and  includes  terms  to 
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allow  for  interfacial  smoothing  and  reindom  “noise”.  The  KPZ  equation  hcis  not  been 
solved  analyticcilly  for  interfaces  growing  in  the  three-dimensional  physical  space,  but 
a  large  body  of  numerical  evidence  shows  that  it  describes  the  evolution  of  a  self-affine 
(fractal)  interface  (see  e.g.,[40]),  with  a  fractal  dimension  that  is  possibly  a  function 
of  varying  growth  conditions  (e.g.,[95]). 

Since  we  aim  to  investigate  how  different  formation  conditions  effect  measurable 
statistical  parameters,  we  construct  a  computer  model  to  simulate  reaction- limited 
kinetics  with  a  tunable  rate  of  growth.  Our  model  is  a  discrete  three-dimensional 
particle-deposition  model  which  is  a  variant  of  the  so  called  “single-step”  model  (SSM) 
proposed  by  [48].  The  SSM  and  variations  of  it  have  been  used  extensively  as  generic 
models  for  interface  growth.  One  reeison  for  this  particular  choice  of  model  is  the 
theoretical  connection  that  can  be  made  with  the  KPZ  equation  [48,  6].  Our  variation 
allows  for  dissolution  to  occur  cis  well  as  deposition,  and  one  can  choose  the  relative 
rates  of  dissolution  versus  deposition  by  changing  the  vcilue  of  a  control  parameter. 
Simulation  results  indicate  that  as  the  ratio  of  dissolution  versus  deposition  at  the 
interface  approaches  unity,  the  fractal  dimension  of  forming  interfaces  increases.  The 
range  of  fractal  dimensions  of  the  simulated  interfaces  lies  between  D  =  2.63  ±  0.005 
and  D  =  2.84  ±  0.01,  in  close  agreement  to  observations. 

We  then  introduce  a  second  variation  of  this  model  in  which  we  allow  the  interface 
to  undergo  peirtial  thermodynamic  equilibration  when  dissolving.  This  allows  for  a 
thermodynamic  distribution  of  dissolution  features,  and  formation  of  etch-pits  and 
holes.  This  second  variation  of  the  model  results  in  non-symmetrical  dissolution  and 
deposition  kinetics,  which  might  prove  to  be  a  more  realistic  description  of  growth 
dyncimics,  since  asymmetrical  functions  for  dissolution  and  precipitation  have  been 
experimentally  observed  for  many  minerals  [65,  53,  54].  Simulations  show  that  statis- 
ticcd  descriptions  of  interfaces  formed  by  Model  II  are  similar  to  those  obtained  from 
Model  I. 

After  studying  how  growth  affects  formation  of  long-range  correlations  on  inter¬ 
faces,  as  measured  via  their  fractal  dimensions,  we  investigate  a  different  geometrical 
property,  the  “roughness  amplitude”,  and  its  relation  to  formation  dynamics.  We  find 
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that  the  amplitude  of  interface  fluctuations  increases  both  with  dissolution  and  the 
fractal  dimension  of  interfaces.  Simulation  results  show  good  qualitative  agreement 
with  our  theoretical  predictions. 

Finally  we  compare  results  from  simulations  to  existing  observations  and  suggest 
an  avenue  for  future  work. 


2.2  Experimental  motivation 

Experimental  measurements,  mostly  motivated  by  a  need  to  better  characterize  and 
predict  properties  of  sedimentary  rocks,  have  shown  that  sedimentary  rocks  are  yet 
another  one  of  the  existing  fractal  objects  to  be  found  in  nature  [5,  30,  94,  84,  21, 
38,  39,  83].  Specifically,  measurements  have  shown  that  pore-grain  interfaces  in  most 
sandstones  and  shales  are  statistically  scale-invariant  over  up  to  four  orders  of  magni¬ 
tude;  from  10~*/im,  the  scale  of  the  smallest  cementing  crystals,  to  the  scale 

of  a  characteristic  pore. 

These  experimental  measurements  were  performed  using  a  variety  of  techniques. 
[21]  covered  the  pore  surfaces  observed  on  thin  sections  with  boxes  of  different  sizes 
to  find  a  power-law  dependence  between  the  size  of  boxes  and  the  number  of  boxes 
needed  to  cover  the  pore  surfaces.  [84]  measured  the  chord-length  distribution  formed 
by  intersections  between  lines  and  pore  surfaces  observed  on  thin  sections  and  frac¬ 
tures.  Autocorrelation  mecisurements  were  also  done  on  thin-sections  [30].  All  the 
above  mentioned  techniques  use  thin-sections  which  are  limited  by  the  polishing  pro¬ 
cess  to  a  resolution  of  l^m.  In  order  to  find  the  statistics  at  the  molecular  level, 
molecular  adsorption  [5]  and  small-angle  scattering  measurements  [94]  were  per¬ 
formed;  these  indicate  fractal  pore  surfaces.  Capillary-pressure  measurements  [18] 
have  been  used  [83]  on  the  intermediate  scale,  between  10"^  and  1/xm  to  provide 
overlapping  data  that  supports  the  continuous  power-law  nature,  over  all  relevant 
scales,  of  pore-greiin  interfaces  in  most  sedimentary  rocks. 

Most  (~  75%)  of  the  fractal  dimension  measurements  of  sedimentary  rock  pore- 
grain  interfaces  presented  by  [38]  (Figure  2-3)  fall  within  the  range  2.6<D<2.8.  [94] 
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found  a  distribution  of  fractal  dimensions  between  2.25  and  2.95.  Figure  2-2  (repro¬ 
duced  from  [85])  shows  thin-sections  from  3  different  sandstones,  with  D  =  2.55, 2.66 
and  2.75.  In  these  thin-sections  one  can  qualitatively  observe  a  trend  of  increasing 
fractal  dimension  with  increasing  amount  of  cementation  due  to  chemical  diagenetic 
alteration  of  pores. 

Figure  2-4  (from  [85])  serves  to  quantify  the  observation  that  D  increases  with 
diagenesis.  When  the  pore-grain  interface  is  fractal,  a  fractal  porosity,  (f>f,  can  be  as¬ 
sociated  with  the  pits  cind  protruding  features  of  these  interfaces.  The  remaining  open 
space,  in  which  one  cein  inflate  an  imaginary  balloon,  is  defined  to  be  the  Euclidean 
porosity,  (j>e  (see  Figure  2-1).  (j>f  is  measured  [85]  using  the  capillary  pressure  method 
of  de  Gennes  [18]  that  predicts  that  the  capillary  pressure  of  a  non-wetting  fluid, 
forced  by  pressure  gradient  to  displace  a  wetting  fluid,  is  indicative  of  the  geometry 
of  pore-grain  interfaces  in  the  rock.  The  total  measured  porosity,  <f)m  =  -f  is  ob¬ 
tained  by  gas-displacement  methods.  The  relative  measure  of  fractal  porosity  versus 
total  porosity,  4>ff4>m  =  ^  constitutes  a  measure  of  diagenetic  alteration.  As 

diagenesis  progresses,  the  total  porosity,  <f>m,  changes  from  being  associated  mciinly 
with  large  open  voids,  <f)e,  to  porosity  associated  mainly  with  pits  and  protrusions  on 
a  rough  interface,  <f)f.  In  Figure  2-4  one  can  see  that  D  increases  with  amount  of 
diagenetic  alteration  as  measured  by 


2.3  Self- affine  interfaces 


In  this  section  we  formulate  a  mathematical  description  of  the  physical  mechanisms 
that  alter  pore-grain  interfaces  in  sedimentary  rocks.  We  study  the  evolution  of  an 
interface  between  two  distinct  regions,  a  region  of  a  pore  filled  with  fluid  and  a  region 
of  solid,  as  shown  in  Figure  2-5.  The  height,  h{x),  of  the  interface  between  the  two 
phases  is  defined  cis  the  distance  of  the  interface  from  a  reference  substrate. 
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2.3.1  Statistical  description 


A  self-affine  interface  is  statistically  similar  to  itself  under  an  affine  transformation; 
i.e.,  when  the  directions  parallel  and  perpendicular  to  the  substrate  are  magnified  by 
different  values  [87,  89]  such  that 

h{x)  «  b-“h{bx),  (2.3) 

where  b  is  the  amount  of  stretching  parallel  to  the  substrate  and  is  the  amount  of 
stretching  perpendicular  to  it.  The  approximation  sign  indicates  that  the  two  sides 
of  the  equation  have  identical  statistical  properties. 

A  useful  statistical  pcirameter  is  the  standard  deviation  of  the  height,  which  will 
be  identified  here,  as  in  much  of  literature  on  interfaces,  as  the  interface  width,  W : 

WiL)  =  {\h{x)-h\y/\  (2.4) 

where  ( )  denote  an  ensemble  average,  and  h  =  1/L^  Jq  h{x)dx  is  the  mean  interface 
height  averaged  over  the  lateral  extent  of  the  system  L.  The  width  of  self-affine 
interfaces  can  be  shown  to  relate  to  the  linear  dimension  of  the  substrate  L  via  a 
power  law  [89]: 

W{L)  ~  i“.  (2.5) 

Here  the  exponent  a  has  a  simple  relation  to  the  fractal  dimension  of  the  interface, 

a  =  3-D,  (2.6) 


for  interfaces  embedded  in  d  =  3. 

Figure  2-6  shows  plots  of  interfaces  (obtained  from  our  simulations  as  outlined 
in  following  sections)  tciken  from  ensembles  with  2  different  fractal  dimensions.  The 
first  one.  Figure  2-6a,  is  a  realization  from  an  ensemble  of  interfaces  with  a  fractal 
dimension  of  D  «  2.63.  It  is  relatively  smooth,  with  greater  dominance  of  long- 
wavelength  structures  over  small-scale  variability.  Figure  2-6b  is  a  realization  taken 
from  an  ensemble  with  D  w  2.84.  This  interface  is  jagged  with  relatively  greater 
short-wavelength  variability. 


21 


After  introducing  the  basic  concepts  used  in  the  study  of  self-affine  interfaces  (for 
a  review  see  [6]),  we  next  examine  continuum  models  of  interfaces  undergoing  generic 
dynamical  processes  of  growth  and  smoothing. 

2.3.2  Continuum  models  for  dynamical  growth  of  interfaces 

Smoothing  processes 

The  time  evolution  of  a  growing  curved  interface  when  subjected  to  smoothing  pro¬ 
cesses  such  as  surface  tension  effects,  surface  diffusion,  mechanical  erosion,  weather¬ 
ing,  recrystalization,  etc,  can  be  described  by  a  simple  diffusion  equation; 

=  „V^h{y.,t)  +  X,  (2.7) 

where  h(x,  t)  is  the  surface  height,  A  is  the  average  velocity  of  interface  growth,  and 
1/  is  an  effective  diffusion  coefficient  from  the  combined  smoothing  effects.  A  diffusion 
term  in  this  context  represents  the  decay  of  a  curved  interface  due  to  processes  that 
favor  a  lower  surface  energy.  Such  processes  only  rearrange  the  solid  phase  mass 
beneath  the  interface  while  conserving  total  mass.  Any  initial  sinusoidal  component 
of  the  interface  which  evolves  by  equation  2.7  decays  exponentially  to  /i(x,  t)  =const 
[13,  p.7-12].  Flat  interfaces,  such  as  those  resulting  from  the  steady-state  solution  of 
the  diffusion  equation  2.7,  correspond  to  a  fractal  dimension  oi  D  =  2  [91]. 

“Noisy  smoothing” 

In  real  rocks  random  fluctuations  can  stem  from  impurities,  anisotropies,  nucleation 
processes,  etc.  We  can  mathematically  model  randomness  in  the  formation  process 
in  a  simple  way  by  letting  the  rate  of  change  in  height  have  uncorrelated  Gaussian 
fluctuations,  T]{x,t),  with  a  zero  mean  and  an  amplitude  of  Q: 

<  7/(x,  t')  >—  2Q6{x  -  x')S{t  -  t').  (2.8) 

Given  such  a  white  noise  structure,  the  full  stochastic  description  of  interfaces 
undergoing  “noisy  smoothing”  is: 

dh 


Interfaces  that  grow  according  to  2.9  have  a  logarithmic  relation  between  the  width 
and  lateral  extent  of  the  system  [20]: 

W\L)  ~  ln(L)  (2.10) 

Hence  an  interface  undergoing  a  “noisy-smoothing”  process  does  not  have  a  power- 
law  dependence  of  W  on  Xt,  as  in  equation  2.5,  and  therefore  does  not  have  a  fractal 
dimension  as  defined  by  equation  2.6.  The  closest  approximation  is  a  fractal  dimension 
of  D  =  d  =  3,  equal  to  the  Euclidean  dimension  in  which  it  is  embedded.  Thus 
equations  2.7  and  2.9  are  insufficient  to  describe  the  dynamical  growth  process  which 
is  responsible  for  formation  of  fractal  interfaces. 

Reaction-limited  growth 

Interface  kinetics  can  also  influence  the  physics  of  growth.  Consider  a  mineral  de¬ 
positing  from  a  saturated  solution  on  a  substrate  such  as  described  in  Figure  2-5. 
The  deposition  process  is  composed  of  two  main  steps:  transport  of  mineral  in  the 
fluid  phase  to  the  solid  interface  via  advection  and  diffusion,  and  incorporation  of 
the  mineral  at  the  interface  via  chemical  reaction.  The  slower  of  the  two  steps  deter¬ 
mines  the  rate  of  interface  growth.  If  reaction  is  slower  than  transport,  the  growth  is 
called  “reaction-controlled”,  while  if  transport  is  slower,  growth  is  termed  “transport- 
controlled”.  Reaction-controlled  kinetics  result  in  disappearance  of  concentration  gra¬ 
dients  in  the  fluid  phase  due  to  the  fast  transport,  so  that  the  concentration  of  mineral 
in  the  fluid  is  constant  in  space  [50].  Thus  the  description  of  growth  simplifies  to  only 
one  variable,  the  height  of  the  interface.  Since  there  are  no  concentration  variations 
in  the  fluid,  every  point  on  the  interface  will  have  equal  probability  to  grow.  Be¬ 
cause  there  are  no  preferred  directions  and  positions  for  growth,  growth  will  occur 
normal  to  the  local  orientation  of  the  interface  and  will  be  statistically  uniform  in 
space.  This  description  is  different  from  transport-limited  growth,  where  one  must 
consider  the  coupling  between  the  interface  height  and  the  mineral  concentration  field 
in  the  fluid  [42].  In  that  case  growth  sites  that  protrude  into  the  fluid  phase  have 
a  higher  probability  for  deposition  than  growth  sites  on  flat  areas,  because  of  the 
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steeper  mineral  concentration  gradient  necir  protrusions.  Figure  2-7  schematically 
illustrates  two  extreme  cases  in  the  evolving  different  shapes  of  interfaces  growing  by 
reaction-controlled  and  diffusion-controlled  mechanisms. 

Recent  experiments  indicate  that  most  sedimentary  rocks  (excluding  some  carbon¬ 
ates  [52,  51])  form  in  the  reaction-limited  regime  [43,  65,  53,  54,  50].  Characteristic 
order-of- magnitude  experimental  reaction  rates  are  JT.  ~  10“®  moles/m^/sec  for  sil¬ 
ica  [65]  and  ~  10“^®  moles/m^/sec  for  kaolinite  [53,  54].  An  order-of-magnitude 
calculation  for  diffusion  rates  of  dissolving  silica  is  [43] 

tr  '^{Ceo  ~ 

Ax)  =  - - - . 

0 

Taking  the  diffusion  rate  constant  V  to  be  10“®  cm^/sec,  the  solubility  of  silica  to  be 
a,  10~®  moles/liter,  the  concentration  at  infinity  to  be  Coo  =  0,  and  the  charac¬ 
teristic  width  of  the  boundary  layer  to  be  5  =  l/xm  ,  then  K-p  ~  10“‘*moles/m^/sec. 
Thus  diffusion  rates  in  such  a  calculation  are  O(IO^)  times  faster  then  reaction  rates, 
resulting  in  a  uniform  concentration  of  mineral  in  the  fluid,  and  reaction-limited 
growth. 

Given  an  interface  that  grows  by  reaction-limited  growth,  (i.e.  normal  to  its  local 
orientation)  with  a  constant  normal  growth  velocity  A,  Figure  2-8  shows,  from  purely 
geometrical  arguments,  that  the  increment  of  growth  projected  onto  the  h  direction 
is 

sh  =  [{xstf  +  (xstvhyy^^ .  (2.11) 

For  small  slopes  this  can  be  written  [29]  as 

^  A  +  iA{VA)^  (2.12) 

where  the  growth  velocity  A  depends  on  the  saturation  of  depositing  minerals  in 
the  fluid  eind  reaction  rates  and  Ccin  be  time  dependent.  Here  we  sheill  cissume  for 
simplicity  that  the  saturation  is  queisi-static,  i.e.  the  saturation  is  effectively  constant 
on  the  time-scales  required  for  the  interface  to  reach  a  statistical  steady-state. 

To  combine  smoothing,  randomness,  and  growth  processes  we  write  an  equation 
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describing  the  evolution  of  an  interface  as 

+  ^X(Vhf  +  ,(x,  t),  (2.13) 

where  we  have  changed  to  a  reference  frame  moving  with  the  average  interface  ve¬ 
locity  A.  Equation  2.13  (commonly  referred  to  as  the  “KPZ”  equation),  was  first 
introduced  for  the  study  of  interface  growth  by  [29].  The  KPZ  equation  has  been 
studied  extensively  as  a  continuum  model  for  evolving  interfaces,  but  due  to  in¬ 
tractable  mathematical  difficulties,  numerical  methods  and  theoretical  investigation 
of  its  statistical  characteristics  constitute  the  main  directions  of  research  [6]. 

The  deterministic  version  of  2.13,  i.e.,  with  =  0,  has  known  solutions.  Re¬ 

sulting  surfaces  develop  as  a  collection  of  paraboloids  joined  together  by  discontinu¬ 
ities  in  Vh.  Normal  growth  results  in  bumps  growing  laterally  as  well  as  sideways,  so 
that  an  interface  with  some  initial  random  configuration  of  features  will  tend  toward 
increasing  dominance  of  long-wavelength  features  over  short-wavelength  features,  as 
seen  in  Figure  2-7a.  The  relaxation  toward  a  flat  interface  in  this  case  is  interesting 
and  quite  different  from  the  ordinary  “diffusion”  dominated  case,  as  described  by  2.7. 
For  example,  in  an  interface  flattening  in  d  =  2,  the  lateral  extent  of  paraboloids  grows 
with  a  power-law  dependence  on  time,  fcister  than  the  decay  of  a  surface  flattening 
due  to  diffusion- like  smoothing  processes  [29]. 

Dynamic  renormalization  group  calculations  predict  that  interfaces  that  evolve 
according  to  the  full  stochastic  equation  2.13  exhibit  statistical  scaling  in  space  and 
time  [29].  The  width  of  these  interfaces  grows  with  time  until  it  reaches  a  steady 
state,  after  which  it  retains  a  constant  value 

VF.(L)  ~  (2.14) 

The  steady-state  interface  is  thus  a  self-affine  fractal,  as  defined  by  2.5  and  2.6.  While 
there  are  no  conclusive  theoretical  predictions,  computer  simulations  in  d  =  3  suggest 
the  existence  of  a  continuous  transition  from  formation  of  non-fractal  interfaces  to 
formation  of  interfaces  with  D  «  2.6  (or  a  0.4)  as  |A|  is  increased  from  0  to  some 
finite  value  [95,  4,  59]. 
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An  intuitive  understanding  of  the  physics  described  by  equation  2.13  can  be  gained 
by  considering  the  effects  of  each  term  in  2.13  on  formation  of  long-range  correlations 
in  the  system: 

a)  When  only  diffusion-like  smoothing  acts  upon  the  interface  {v  ^0  eind  A,  77  =  0  in 
2.13,  resulting  in  2.7),  the  steady-state  interface  is  h{x.)  =const,  and  since  all  points 
have  identical  height,  a  “height-height”  correlation  function  does  not  decay  in  space 
and  D  =  2. 

b)  When  white  noise  is  added  (i/,?;  ^  0,  A  =  0,  resulting  in  2.9),  the  never-decaying 
correlations  obtained  in  case  (a)  are  diminished  by  the  noisy  reindom  forcing.  This 
results  in  an  interface  width  W  that  grows  only  logarithmically  with  system  size 
corresponding  to  the  limit  a  0  (D  — »  3)  in  equation  2.6. 

c)  When  reaction-limited  growth  is  also  present  ^  0,  resulting  in  the  full 

equation  2.13),  “bumps”  grow  normal  to  the  interface;  thus  when  they  grow  upwards 
they  also  grow  sideways  (Figure  2-7a)  at  a  rate  feister  than  diffusion  would  predict, 
allowing  local  “height  information”  to  be  transmitted  laterally.  Hence  normal  growth 
enhances  formation  of  long-range  correlations  and  long-wavelength  features,  while 
suppressing  or  smoothing  out  short  wavelength  features,  and  thus  decreases  the  fractal 
dimension  of  the  forming  interface. 

Equation  2.13  therefore  represents  a  balance  between  factors  (diffusion  and  re¬ 
action  limited  growth)  which  tend  to  reduce  small-scale  features  (decrease  fractal 
dimension)  cind  a  factor  (reindom  forcing)  which  tends  to  relatively  increase  small- 
scale  features.  The  balance  that  is  struck  by  the  coefficients  of  2.13  should  determine 
the  fractal  dimension  of  the  pore  interfaces. 

A  common  method  for  quantitative  study  of  growing  interfaces  is  the  construction 
of  simple  discrete  models  governed  by  processes  similaj  to  those  described  by  the  re¬ 
spective  continuous  equations.  This  approach  avoids  the  severe  sensitivity  that  direct 
numerical  solutions  to  2.13  exhibit.  Such  analog  discrete  models  are  also  appealing 
due  to  their  relatively  simple  implementation  and  the  fact  that  in  most  interface 
growth  problems  (eis  in  our  problem  of  growth  of  crystals  on  interfaces)  the  physical 
system  studied  is  actually  discrete  by  nature.  Nevertheless,  continuous  equations  such 
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as  2.13  provide  a  predictive  physical  framework  for  the  discrete  studies.  Here,  the 
continuous  representation  is  useful  both  for  isolating  the  different  processes  involved 
in  formation  of  pore-grain  interfaces  and  for  providing  some  physical  insight  into  dy¬ 
namics  of  formation  of  long-range  correlations.  For  a  review  of  recent  approaches  and 
results  in  the  study  of  fractal  interface  growth  see  [45]  and  [6]. 


2.4  Computer  simulations 

We  next  present  two  simple  discrete-particle  models  of  interfaces  roughening  by  de¬ 
position  and  dissolution,  variations  of  the  so-called  ’single-step’  model  (SSM)  [48]  . 
The  average  properties  of  the  SSM  can  be  calculated  and  shown  to  correspond,  to  a 
first  approximation,  to  equation  2.13  [48,  60].  It  is  this  theoretical  correspondence,  as 
well  as  the  existence  of  a  physical  analog  between  mechanisms  in  rocks  and  deposition 
and  dissolution  in  the  model,  that  led  us  to  choose  the  SSM  over  the  multitude  of 
other  discrete-particle  models  used  to  study  the  KPZ  equation. 

The  original  SSM  (Figure  2-9)  starts  with  a  square  lattice  that  is  filled  by  steps 
in  a  checker-board  manner;  i.e.,  every  filled  site  is  a  step  of  height  1  surrounded  by 
nearest  neighbor  holes  of  height  0.  At  each  successive  time  step,  a  site  is  chosen  at 
random  from  all  the  sites  that  are  local  minima  (i.e.,  sites  that  are  lower  than  any  of 
their  nearest  neighbors).  The  chosen  site  is  then  filled  by  a  block  of  height  2,  so  that 
it  now  becomes  a  local  maximum,  1  step  higher  than  its  neighbors.  Qualitatively,  the 
SSM  captures  the  three  generic  physical  processes  described  by  the  KPZ  equation  in 
the  following  way: 

I)  The  SSM  has  an  intrinsic  smoothing  process,  corresponding  to  a  diffusional  term 
V  in  2.13,  due  to  the  requirement  that  the  choice  of  deposition  sites  must  be  among 
the  ones  that  are  local  minima,  thus  effectively  “smoothing”  away  holes. 

II)  Randomness,  corresponding  to  77  in  2.13,  is  incorporated  in  the  SSM  by  the  ran¬ 
dom  choice  among  all  available  sites. 

III)  Reaction-limited  growth  is  incorporated  because  growth  is  restricted  by  the  avail¬ 
ability  of  growth  sites,  rather  than  by  the  supply  of  blocks  from,  say,  a  diffusing  field. 


27 


The  non-linear  term  in  the  continuum  description  of  the  SSM  emerges  from  a  geo¬ 
metrical  argument  similar  to  the  one  made  in  deriving  2.13.  In  Figure  2-8  the  change 
in  local  height  (5A)  during  normal  growth  is  an  increasing  function  of  the  local  slope 
(VA),  so  that  A  >  0  in  2.13.  In  a  2d  SSM,  a  locally  flat  interface  has  a  height  config¬ 
uration  of  h{xi)  =  c  for  2  even,  h(xi)  =  c  -|- 1  for  i  odd,  where  c  is  a  constant.  In  this 
case,  half  the  sites  are  local  minima  and  are  available  for  growth,  and  the  interface 
cam  grow  rapidly.  On  the  other  hand,  for  an  inclined  interface  {h(xi)  =  a  +  bi)  there 
axe  no  sites  available  for  growth,  since  no  site  is  a  local  minimum.  Hence,  as  with 
the  KPZ  equation,  the  SSM  model  enforces  a  dependence  of  local  height  ch2Lnge  (Sh) 
on  local  slope,  but  in  the  opposite  sense  (i.e.  A  <  0)..  A  quantitative  derivation  of 
2.13  from  the  average  properties  of  the  SSM  was  made  using  mappings  of  the  SSM 
to  Ising  spin  and  lattice  gas  models,  and  can  be  found  in  [48]  and  [6]. 

2.4.1  Model  I:  Symmetric  dissolution  and  precipitation 

Specifications  We  are  interested  in  testing  the  hypothesis  that  the  fractal  dimen¬ 
sion  of  forming  interfaces  depends  on  the  relative  amplitudes  of  noise  tj,  smoothing 
rate  t/,  and  reaction- controlled  growth  rate  A.  We  propose  to  control  the  reaction  con¬ 
trolled  growth  rate  by  allowing  dissolution  of  blocks  in  dynamics  that  mirror  those  of 
deposition.  At  each  time-step  a  deposition  event,  as  described  above,  will  occur  with 
a  probability  (0  >  p+  >  1)  and  a  dissolution  event  will  occur  with  a  probability 
p_  =  1  —  p+.  A  dissolution  event  is  defined  to  be  the  subtraction  of  a  block  of  length 
2  from  a  site  randomly  chosen  cimong  edl  the  sites  that  are  local  maxima.  By  allow¬ 
ing  particles  to  attach  and  detach  to  the  interface,  we  hope  to  simulate  molecular 
exchange  across  phase  boundaries;  increasing  p_  increases  the  number  of  particles 
leaving  the  interface  versus  the  number  attaching  to  it. 

Two-dimensional  analogs  of  our  model  have  been  studied  by  [60]  and  [6].  They 
calculate  that,  on  average,  the  evolution  of  the  simulated  interface  is  described  by 
2.13.  Parameters  v  and  Q  cure  constants  independent  of  p+,  while 

=  -(P+  -  P-)-  (2.15) 
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Qualitatively,  2.15  describes  a  relation  between  the  average  growth  velocity  of  the 
interface  (oc  (p+  —  p_))  and  the  non-linear  coefficient  A.  As  expected  from  analog 
calculations  for  the  continuum  model  (Figure  2-8  and  equation  2.12)  the  magnitude 
of  the  non-linear  coefficient  increases  with  increasing  interface  growth  velocity.  Here 
|A1  is  maximum  when  p+  =  1,  and  decreases  to  0  when  dissolution  balances  deposition 
and  the  interface  has  no  net  growth. 


Results.  Simulations  of  growing  interfaces  were  performed  for  different  system  sizes , 
with  p+  varying  from  0.5  to  1.  (Note  that  p+  =  0.5  is  a  symmetry  point,  and  results 
obtained  for  adveincing  interfaces  are  applicable  to  retreating  interfaces,  with  p_  and 
p+  exchcinged).  Representative  interfaces  with  p+  =  1  and  p+  =  0.6  are  shown  in 
Figure  2-6.  To  quantitatively  test  whether  resulting  interfaces  are  self-affine,  the 
interface  width,  W{t)  (as  defined  by  equation  2.4)  is  measured  for  different  system 
sizes  L  and  averaged  over  an  ensemble  of  300  simulations  performed  with  different 
random  numbers.  We  find  that  the  width  of  interfaces  grows  as  a  function  of  time 
until  it  reaches  a  statistically  constant  saturation  value,  W,.  Thereafter  it  exhibits 
a  power-law  dependence  on  system  size,  L,  and  remains  self-affine,  obeying  equation 
2.14.  We  ascribe  this  behavior  to  growth  of  ‘bumps’  both  vertically  and  horizontally, 
as  explained  for  Figure  2-7a,  until  a  saturation  value  for  the  amplitude  of  the  largest 
bumps  is  obtained  when  the  wavelength  corresponding  of  the  lateral  extent  of  the 
largest  features  reaches  the  system  size.  The  initial  transient  phase  of  power-law 
growth  in  time  is  chciracteristic  to  interfaces  that  obey  dynamics  described  by  2.13. 
The  initial  power-law  growth  of  our  model  agrees  well  with  theoretical  predictions 
[40]  for  all  p+^O.d.  Although  this  transient  evolution  of  interfaces  is  an  interesting 
aspect  of  the  problem,  we  limit  the  discussion  in  this  paper  to  the  non-equilibrium 
steady-state,  since  that  is  where  we  can  make  comparison  with  experiments. 

To  demonstrate  the  fractal  nature  of  the  resulting  statistically  steady-state  inter¬ 
faces,  we  plot  logio  W,  versus  log^  L.  For  self-affine  interfaces  equations  2.14  and  2.6 
should  hold,  and  hence  we  expect  that  the  width  of  self-affine  interfaces  will  plot  as 
straight  lines  on  this  graph  with  a  slope  which  is  equal  to  3  -  D.  Figure  8a  shows 
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such  plots  for  various  values.  Note  that  cis  p+  decreases,  the  slope  of  the  graph 
decreases,  indicating  an  increase  in  the  fractal  dimension  of  the  resulting  interfaces 
with  increasing  amovmt  of  dissolution.  For  p+;$0.6  the  results  of  simulations  are 
not  well-fit  by  a  straight  line;  these  interfaces  are  not  self-affine.  Figure  2-lOb,  a 
plot  of  versus  Ini,  demonstrates  that  2.10  provides  a  better  fit  for  the  data  for 
high  dissolution  rates.  This  is  consistent  with  the  expectation  that  a  transition  from 
power-law  to  ~  Ini  behavior  will  occur  when  p_  — >  p+  and  A  ^  0. 

After  obtaining  D  for  various  p+  values  from  Figure  2-lOa,  we  investigated  func¬ 
tional  relations  between  deposition  rate,  dissolution  rates,  and  fractal  dimensions  of 
forming  interfaces.  In  Figure  2-11  we  plot  D  versus  a  nondimensional  pcirameter 
which  we  term  the  “dissolution-deposition  ratio”.  The  dissolution-deposition 
ratio  measures  the  rate  of  pairticles  leaving  the  interface  (2p_)  versus  the  rate  of 
particles  attaching  to  it  (2p+). 

The  fractcil  dimension  of  resulting  interfaces  is  seen  to  increase  with  the  dissolution- 
deposition  ratio,  with  a  curve  showing  that  for  p_/p+  =  0,  £)  =  2.63  ±  0.005,  in 
agreement  with  previous  simulations  of  the  pure  deposition  SSM  [48],  and  as  p_/p+ 
increeises,  D  approaches  a  limiting  value  of  2.84±0.01.  An  increase  in  the  dissolution- 
deposition  ratio  above  p_/p+  ~  0.7  results  in  suspected  loss  of  fractal  behavior  and 
a  transition  to  logarithmic,  rather  than  power-law,  dependence  of  W,  on  L. 


2.4.2  Model  II:  Asymmetric  model  of  dissolution  and  pre¬ 
cipitation 

Specifications.  We  next  construct  a  Vciriation  of  Model  I  that  models  dissolution 
in  partial  thermodynamic  equilibrium  [93].  This  is  done  in  order  to  investigate  the 
effects  of  a  finite  probability  for  the  development  of  features  such  as  etch  pits  and 
holes  that  can  be  found  in  rocks.  Allowing  for  partial  thermodynamic  equilibration 
in  the  dissolution  step,  but  not  in  the  deposition  step,  creates  an  asymmetry  between 
the  two  processes.  Asymmetrical  functions  for  dissolution  and  precipitation  have  been 
experimentally  observed  for  many  minerals  [65,  53,  54]. 
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In  Model  II,  deposition  occurs  with  probability  p+  identical  to  the  deposition  step 
in  Model  I,  but  a  temperature-dependent  dissolution  step  occurs  with  probability 
P_  =  1  _  in  a  way  that  allows  for  thermodynamic  equilibration  of  the  interface. 
The  dissolution  step  is  constructed  as  following:  A  site  i  is  randomly  chosen  among 
all  the  sites  of  the  interface  and  the  local  surface  area 


^  6 


is  measured,  where  summation  is  over  the  4  nearest  neighbors.  If  a  block  is  to  be  dis¬ 
solved  at  site  t,  the  change  in  surface  area,  AEi  =  -  2  -  hj+il  —  \hi  -  hi+s\), 

is  ccilculated.  The  probability  q  of  occurrence  of  a  dissolution  event  at  site  i  is  then 
defined  to  be 


9  = 


1 

g-AElkT 


if  AJ5  <  0 
if  >  0 


(2.16) 


Here  kT  is  a  tunable  model  temperature,  roughly  analogous  to  a  thermodynamic 
temperature.  Dissolution  at  site  i  will  thus  happen  if  a  dissolution  event  results  in 
reduced  or  equivalent  surface  area,  and  will  have  an  exponentially  decaying  probability 
to  occur  if  the  surface  area  is  increased  by  the  dissolution  step.  Finite  probability  for 
increcising  the  surface  2U'ea  is  allowed  in  order  to  model  etch- pits  and  holes,  features 
observed  in  rocks.  If  dissolution  of  a  block  did  not  occur  at  the  site  i  first  chosen,  a 
different  site  is  randomly  chosen  eunong  all  the  sites  of  the  interface  and  a  dissolution 
event  is  attempted  (according  to  rule  2.16)  at  the  newly  chosen  site,  and  so  on,  until 
an  attempt  to  dissolve  is  successful.  The  only  effect  of  these  repetitive  attempts  at 
dissolution  is  to  force  p_  to  be  constant  and  equal  to  1  —  p+.  At  the  zero  temperature 
limit  this  model  reduces  to  Model  I.  At  high  enough  temperatures  the  dissolution 
step  introduces  only  noise  to  the  system  while  reducing  the  growth  rate. 

Because  dissolution  introduces  a  thermodynamical  equilibration  procedure  that 
is  not  duplicated  in  the  deposition  event,  the  dyncimics  of  retreating  interfaces  with 

<  0.5  in  Model  II  are  not  the  mirror  image  of  interfaces  with  p+  >  0.5.  This 
asymmetry  is  the  fundamental  difference  between  Model  II  and  model  I. 
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Results.  We  performed  statistical  studies  for  interfaces  forming  at  dilferent  kT,p+, 
and  L.  At  all  temperatures  studied  we  find  that,  similar  to  the  results  of  Model 
I,  the  width  of  interfaces  formed  by  Model  II  have  a  transient  stage  of  dynamical 
scaling,  after  which  they  reach  a  statisticed  steady-state.  Graphs  of  \ogW,{L)  as 
function  of  logL  eure  plotted  in  Figure  2-12a  and  b  for  kT  =  1  and  kT  =  100, 
respectively.  It  is  demonstrated  that  power-law  models  provide  a  good  fit  to  the 
data,  although  logarithmic  dependence  between  the  width  and  size  of  the  system  for 
high  dissolution  rates  at  low  temperatures  cannot  be  ruled  out.  The  fractal  dimension 
for  kT  =  1  (Figure  2-12a)  is  seen  to  increase  with  incretisingp_/p+  similarly  to  results 
of  Model  I.  For  kT  =  100  (Figure  2-12b)  we  note  a  different  behavior  them  for  the  low 
temperatures.  The  fractail  dimension  (deduced  from  the  measured  slope)  of  interfaces 
is  nearly  constant  with  increased  dissolution,  but  the  amplitude  of  the  width  increases 
with  decreasing  p+. 

Figure  2-13  shows  the  fractal  dimension  of  interfaces  as  function  of  p_/p+.  The 
four  curves  correspond  to  four  different  temperatures:  kT  =  10“*,  1,10*  from  this 
model,  and  kT  =  0  replotted  from  Model  I.  At  high  temperatures  {kT  =  10*)  the 
fractal  dimension  is  nearly  constant  as  a  function  of  relative  dissolution  rate.  For 
low  temperatures  all  curves  follow  the  same  trend  of  increased  D  with  increased 
dissolution.  We  believe  that  interfaces  formed  at  low  temperature  are  not  fractal  for 
p_/p+>0.7. 

2.4.3  Roughness  of  interfaces 

The  notion  of  roughness  of  an  interface,  i.e.  the  eimplitude  of  its  fluctuations,  cam  be 
qucmtified  by  measuring  the  prefactor  in  equation  2.14,  now  rewritten  as 

W,  =  A{D)L^-^  (2.17) 

The  prefactor,  A(D),  is  termed  the  “roughness  amplitude”,  amd  measures  the  ampli¬ 
tude  of  the  undulations  of  interfaces. 

In  Figure  2-14  we  plot  A,  Ccilculated  for  all  temperatures,  from  intercepts  of  lines 
in  Figures  2-lOa  and  2-12a,b  with  the  width  axis,  versus  p_/p+.  We  find  that  the 
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roughness  amplitude  increcises  approximately  linearly  with  the  dissolution-deposition 
ratio.  It  is  interesting  to  note  that  data  from  all  temperatures  follow  the  same  lin¬ 
ear  trend,  except  where  interfaces  are  not  likely  to  be  fractal.  At  low  temperatures, 
where  dissolution-deposition  ratio  is  high  enough  (p_/p+^0.7),  the  roughness  reaches 
a  plateau  eind  diverges  from  the  linear  trend.  We  attribute  this  to  the  logarithmic 
dependence  of  on  L  for  dissolution-deposition  ratios  close  to  iinity.  Thus  we  can 
use  the  divergence  from  the  linear  trend  in  Figure  2-14  as  another  indication  for  cir¬ 
cumstances  for  which  interfaces  cannot  be  well  modeled  eis  fractcds.  The  dissolution- 
deposition  ratio  obtedned  for  this  divergence  (p_/p+>0.7)  is  consistent  with  that 
determined  from  Figures  2-10  and  2-12. 

The  dependence  of  A  on  D  can  be  predicted  for  the  zero  temperature  model  by 
using  two  constraints.  The  first  constraint  emerges  from  the  “single-steppedness”  of 
Model  I:  the  square  of  the  slope  of  the  height  at  any  given  site  must  always  be  equal  to 
unity  (j/it  —  hj+i P  =  1).  The  second  constraint  is  that  the  power  spectra  of  self-affine 
interfaces  have  a  power-law  form  [87).  The  details  of  predicting  A{D)  are  given  in 
Appendix  A. 

Figure  2-15  shows  that  for  fcT  =  0,  A  is  an  increasing  function  of  D,  which  means 
that  on  short  length-sccdes  interfaces  with  high  fractal  dimensions  appear  rougher 
and  have  larger  undulations  than  interfaces  with  lower  fractal  dimensions.  The  dis¬ 
crepancy  between  simulation  results  and  theoretical  predictions  is  probably  due  to 
the  fact  that  the  power  spectrum  in  the  theoretical  calculation  of  A{D)  is  assumed  to 
follow  a  power-law  at  all  wavelengths  (as  given  in  equation  A. 9),  but  in  reality  only 
follows  this  behavior  between  high  and  low  wavenumber  cutoffs. 

We  note  that  although  A  is  an  increcising  function  of  p_/p+  and  D  (Figures  2-14 
and  2-15),  the  measured  width  W  of  interfaces  of  large  enough  lateral  extant  A  is  a 
decretising  function  of  p+lp-  and  D.  This  is  because  as  the  system  size  increases, 
W  increases  as  well,  but  more  slowly  for  interfaces  with  high  fractal  dimensions  than 
for  interfaces  with  low  fractal  dimensions,  which  cam  be  seen  from  both  equation  2.17 
eind  Figure  2-10.  Thus,  for  a  “system-sized  elephant”  an  interface  with  a  high  fractal 
dimension  appears  smoother  than  one  with  a  lower  fractal  dimension,  while  for  a 


33 


“particle-sized  ant”  an  interface  with  a  high  fractal  dimension  is  rougher  than  one 
with  a  low  fractcd  dimension. 


2.5  Summary  and  discussion 

2.5.1  Summary  of  model  and  results. 

Motivated  by  experimental  data  indicating  fractal  pore  surfaces  in  sedimentary  rocks, 
we  have  developed  an  analytical  description  and  a  simple  computer  model  for  reaction- 
controlled  growth  of  interfaces.  Analytical  arguments  lead  to  the  KPZ  equation 
2.13  [29],  a  non-lineair  partial  differentieil  equation  extensively  studied  as  a  model 
for  growth  of  self-affine  interfaces  (e.g.,  [6]).  This  equation  describes  dynamical  in¬ 
terface  evolution  governed  by  diffusion-like  smoothing,  reaction-limited  growth,  and 
random  events.  Our  goal  in  studying  such  dyneimical  descriptions  of  interface  growth 
was  to  find  a  link  between  physical  processes  that  govern  growth  and  geometrical 
properties  of  resulting  interfaces.  Such  a  link  can  help  constrain  formation  history  of 
rocks  by  measuring  their  geometrical  properties. 

In  order  to  study  how  different  dynamical  processes  affect  the  steady-state  statis¬ 
tics  of  interfaces  we  have  constructed  a  discrete  particle  deposition  and  dissolution 
model  which  incorporates  reaction-limited  growth,  interfacial  smoothing,  and  random 
“noise”  processes  at  a  growing  interface  and  provides  a  control  over  the  relative  rates 
of  these  processes.  The  average  properties  of  interfaces  formed  by  this  model  were 
shown  (in  d  =  2)  to  correspond  to  the  KPZ  equation.  We  qualitatively  explain  this 
correspondence. 

Interfaces  formed  by  our  model  go  through  a  transient  stage  of  roughening  after 
which  they  reach  a  statistical  steady  state,  where  the  interfaces  still  grow,  but  their 
statistical  characteristics  remaiin  constant.  The  steady-state  interfaces  are  fractcil  for 
most  parameter  ranges,  with  the  fractal  dimension  increasing  from  2.63  ±  0.005  to 
2.84  ±  0.01  as  the  dissolution-deposition  ratio  is  increeised.  For  dissolution-deposition 
ratios  approaching  unity,  interfaces  formed  are  no  longer  fractal.  These  results  are 
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consistent  with  am  expected  transition  from  fractal  to  non-fractal  interfaces  when  the 
magnitude  of  the  non-lineeir  term  in  the  KPZ  equation  is  decreased  from  a  finite  value 
to  0.  The  nature  of  the  transition  is  not  theoretically  predicted. 

We  also  find  that  the  “roughness  amplitude”  of  interfaces,  A  (as  defined  in  equa¬ 
tion  2.17),  increcises  with  increasing  dissolution-deposition  ratio  and  increasing  fractal 
dimension  of  simulated  surfaces.  This  behavior  is  in  relatively  good  agreement  with 
our  theoreticcil  predictions. 

2.5.2  Comparison  with  experiments 

The  emerging  physical  picture.  Laboratory  measurements  show  that  the  ma¬ 
jority  of  sedimentary  rocks  are  fractal  with  2.55;SJ9;$2.8,  as  seen  in  Figures  2-3  and 
2-4.  This  range  coincides  approximately  with  our  simulation  results.  The  experi¬ 
mental  observations  also  show  that  rocks  with  highly  diagenetically  altered  porosity 
(generally  the  samples  with  higher  content  of  cementing  materials,  and  more  evi¬ 
dence  of  dissolution  and  precipitation)  correspond  to  the  higher  fractal  dimensions. 
We  propose  the  following  scenario  to  expleiin  the  observed  trend  in  the  geometry 
of  pore-grain  interfaces:  Necir-surface  sedimentary  rocks  are  usually  part  of  a  large 
scale  system  through  which  fluid  is  flowing  at  non-negligible  flow  rates.  These  rocks 
are  thus  prevented  from  reaching  global  chemiceil  equilibrium  as  one  would  expect 
for  samples  in  a  closed  system.  Diagenetic  processes  occur  in  the  forming  rock  as 
a  consequence  of  this  non-equilibrium  situation.  Since  diagenesis  generally  acts  to 
reduce  permeability  [69,  63],  the  rock  becomes  more  resistant  to  fluid  flow  with  time. 
Although  global  equilibrium  is  not  reached,  the  reduction  in  flow  rates  results  in  pore 
fluids  spending  more  time  in  a  pore  and  thus  becoming  more  locally  chemically  equi¬ 
librated  with  surrounding  solid.  Thus,  growth  and  dissolution  ionic  fluxes  at  the  pore 
surface  begin  to  equilibrate  cind  p_/p+  increases.  Figure  2-11  shows  that  the  fractal 
dimension  increaises  when  p-fp+  increases,  while  Figures  2-2  eind  2-4  show  that  the 
fractal  dimension  incre^es  with  diagenetic  alteration.  The  model  and  experiments 
together  suggest  that  growth  and  dissolution  in  a  finite  volume  lead  to  a  unique  “di¬ 
agenetic  pathway”  that  is  descriptive  of  pore  evolution  in  many  similar  rocks.  At 
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the  beginning  of  the  pathway  the  pore  is  relatively  open,  crystal  growth  is  reaction 
limited  with  a  small  p_/p+  (or  small  p+/p_,  in  the  symmetrical  case  of  dissolving 
interfaces)  and  the  resulting  fractal  dimension  is  near  2.6.  As  the  pore  space  is  filled, 
the  balance  of  growth  and  dissolution  rates  shifts  toward  unity,  p_/p+  increases,  cind 
so  the  frcictal  dimension  increases.  Following  this  model,  [85]  find  a  limiting  fractal 
dimension,  D  =  2.75,  reached  by  the  competition  between  rates  of  ionic  diffusion  on 
an  increasingly  rough  surface,  and  reaction. 

Finally,  it  is  possible  to  tie  the  physical  picture  with  a  recent  resxilt  that  demon¬ 
strates  (Figure  2-16)  that  the  permeability,  /C,  is  related  to  and  /«,  the  Euclidean 
porosity  and  length  scale,  and  not  to  the  total  porosity  and  length  scales  [85,  2],  via: 

x:  =  (2.18) 

Since  <f>f/^m.  =  1  —  <l>e/<t>m.  increases  with  diagenesis  (see  Figure  2-4),  equation  2.18 
predicts  that  permeability  will  decrease  with  time,  even  if  <f)m  stayed  constant  (a 
prediction  in  agreement  with  experimental  results  in  [69]).  The  reduced  fluid  flow 
rates  result  in  more  equilibration  and  so  the  fractal  dimension  will  increase,  and 
(f)f  and  (f)e  will  consequently  change.  Thus,  the  evolving  microscopic  features  on 
pore-grain  interfaces  both  influence  and  are  influenced  by  the  large  scale  transport 
properties. 

Compeirison  of  model  results  for  the  “roughness  amplitude”  with  data  from  real 
rocks  should  be  approached  with  caution  because  of  the  large  number  of  variable 
parameters.  These  can  be  dealt  with  by  constructing  suitable  transformations  (e.g. 
doubling  molecular  size  would  result  in  double  roughness  amplitude).  Although  at 
this  point  we  have  no  reliable  data  to  which  we  ctin  compare  our  roughness  results, 
estimates  of  A  from  [38]  show  a  trend  of  increasing  A  with  D,  cis  predicted  quan¬ 
titatively  from  our  models.  Our  simulations  also  predict  that  diagenetically  altered 
pores  that  appecir  rough  on  the  crystal  scale  will  appear  smooth  on  the  pore  scale, 
while  rocks  which  have  low  D  and  little  diagenetic  alteration  will  appear  smooth  on 
the  crystal  scale  but  rough  and  strongly  undulating  on  the  pore  scale.  At  this  point 
we  do  not  have  enough  data  to  check  this  prediction. 
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Deviations  of  some  observations  from  the  predicted.  Simulation  results  do 
not  predict  the  observations  of  sedimentary  rocks  with  D^2.6.  Since  it  is  our  in¬ 
tention  to  capture  the  dominant  physical  processes  in  pore-greiin  interfaces,  why  can 
we  not  simulate  these  existing,  though  less  common,  observations?  Most  likely  our 
model  does  not  adequately  describe  till  natural  growth  environments,  in  particular 
the  different  growth  mechanisms.  We  propose  that  the  small  percentage  of  rocks  that 
have  a  fractal  dimension  which  cannot  be  explained  by  our  model  should  serve  as  a 
test  for  a  point  of  departure  of  the  formation  conditions  from  the  ones  assumed  by 
our  models.  For  example,  transport-controlled  growth,  which  was  not  investigated 
here,  may  produce  completely  different  results  from  our  model.  Transport-limited 
growth  might  be  the  mode  of  growth  for  some  carbonates  (e.g.,  calcite  and  arogonite 
at  certain  pH  levels  [52,  51])  that  are  highly  soluble.  Bedford  limestone,  a  carbonate 
for  which  D  =  2.35  [39],  serves  as  an  example  for  transport  limited  growth  leading 
to  fractal  interfaces,  with  a  fractal  dimension  quite  different  than  that  predicted  by 
our  model.  Another  possibility  is  that  all  the  processes  that  we  have  termed  “noise” 
are  not  uncorrelated  as  we  postulate.  While  uncorrelated  noise  might  be  a  good 
assumption  in  most  Ccises  for  forming  sedimentary  rocks  (due  to  the  generally  short 
range  nature  of  the  forces  exerted  by  ions  on  the  interface,  the  random  position  of 
impurities  and  orientation  of  grains  on  which  growth  occurs,  the  random  process  of 
nucleation,  and  a  variety  of  other  conditions),  one  can  imagine  cases  where  random 
events  tend  to  be  correlated  in  space  and  time,  such  as  when  one  mineral  acts  to  lower 
the  surface  energy  for  a  second  mineral  to  crystallize,  or  when  events  are  correlated 
in  space  by  certain  directions  of  growth  being  energetically  preferred.  By  introducing 
power-law  correlated  noise,  interfaces  may  be  formed  with  a  continuously  varying 
fractal  dimension  between  2  and  3,  as  demonstrated  by  computer  and  theoretical 
models  [49,  46,  47].  As  one  might  naively  expect,  forcing  external  correlations  (anti¬ 
correlations)  on  the  growth  process  indeed  increases  (decrecises)  correlations  between 
points  on  the  interface  and  results  in  a  lower  (higher)  surface  fractal  dimension. 
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2.6  Conclusion 


Our  proposed  physical  and  numericed  model  addresses  the  three  requirements  posed 
in  the  introduction:  it  is  a  model  of  non-equilibrium  growth,  it  is  independent  of 
mineralogy,  and  it  agrees  with  the  observations  that  higher  fractal  dimensions  are 
found  in  rocks  that  are  more  diagenetically  altered.  Thus,  we  believe  that  the  general 
physical  mechanism  of  growth  in  most  shales  and  sandstones  can  be  captured  by  the 
simple  processes  of  reaction-limited  growth  eind  smoothing  in  a  noisy  system.  This 
microscopic  growth  is  in  turn  linked  to  the  large  scale  permeability  changes  occurring 
during  diagenesis. 

This  work  constitutes  a  first  step  in  using  geometrical  constraints  to  study  the 
dynamical  history  of  the  formation  of  rocks.  More  quantitative  observations  of  geo- 
metriccil  properties  as  well  as  more  experiments  for  controlled  growth  in  the  laboratory 
are  necessary  before  a  comprehensive  theory  can  be  developed. 


38 


Figure  2-1:  Schematic  diagram  of  a  single  pore  in  a  sedimentary  rock.  Pore-grain 
interfaces  in  sedimentary  rocks  are  genereJly  quite  convoluted  with  geometrical  struc¬ 
tures  formed  by  cementing  crystcds  and  corroded  etch  pits  and  holes.  The  total 
porosity,  <l>m,  is  a  sum  of  ^e>  the  porosity  associated  with  Euclidean  open  pore  space, 
and  (f)f,  the  porosity  associated  with  fractal  undulations  of  the  pore-grain  interface. 


Figure  2-2:  From  [85].  Thin  sections  of  three  seindstones.  (A)  Table  sandstone  with 
a  fractal  dimension  of  D  =  2.55.  (B)  Price  River  sandstone  with  D  =  2.66.  (C) 
Coconino  sandstone  with  D  =  2.75.  The  fractal  dimension  increases  with  increcising 
volume  of  cementing  material  and  dissolution  of  the  initial  sand  grains. 
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Figure  2-3:  A  histogram  of  fractcil  dimensions  for  27  different  rocks.  Bins  axe  of  equal 
size.  The  data  are  from  [38].  75%  of  the  observed  rocks  have  a  fractal  dimension 
between  2.6  and  2.8. 
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Figure  2-4:  From  [85].  The  measured  fractal  dimension  D  of  pore-grain  interfaces 
versus  the  porosity  cissociated  with  fractal  surfaces  over  the  measured  total  porosity, 
<f>fl4>Tni  of  11  different  sandstones  and  1  shale.  The  plot  shows  a  monotonic  increzise 
of  D  with 


42 


Figure  2-5:  A  schematic  view  of  the  interface  between  a  fluid-filled  rock  pore  and  a 
solid  rock  material.  This  is  the  setting  of  the  theoretical  problem:  an  interface,  with 
no  overhangs,  constitutes  the  boundary  between  fluid  and  solid  phases.  The  interface 
position  is  described  as  a  height  deviation,  h{x,t),  from  a  substrate.  L  is  the  lateral 
extent  of  the  substrate  on  which  the  interface  grows. 
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Figure  2-6:  Plots  of  two  interfaces  tciken  out  of  ensembles  with  2  different  fractal 
dimensions.  The  interfaces  were  obtained  from  simulations  of  Model  I,  performed 
on  lattices  of  size  96  x  96.  Deirk  shadings  indicate  lower  then  average  height  and 
light  colors  indicate  heights  greater  then  average,  (a)  A  sample  realization  from 
am  ensemble  with  a  fractal  dimension  of  =  2.63.  It  shows  relatively  subdued 
short- wavelength  features  (formed  with  =  1).  (b)  A  sample  realization  from  an 
ensemble  with  D  =  2.84.  This  interface  is  jagged  with  relatively  more  power  to  short 
wavelength  features  (formed  with  p+  =  0.6). 
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Figure  2-7:  A  schematic  drawing  of  growth  by  (a)  only  reaction  controlled  kinetics 
and  (b)  only  transport  controlled  kinetics  (shown  in  the  limit  of  infinitely  fast  reaction 
rates).  Both  interfaces  start  from  the  same  initicil  condition,  and  successive  profiles 
correspond  to  propagation  in  time,  (a)  is  growing  normal  to  its  loccil  orientation, 
with  a  constant  normal  growth  rate.  Large  bumps  grow  at  the  expanse  of  small  ones, 
creating  parabolas  which  are  joined  by  discontinuities  in  Vh  [29].  (b)  illustrates  the 
Mullins- Sekerka  instability.  Protruding  features  create  steep  concentration  gradients 
in  the  fluid  phase,  thus  increasing  transport  of  mineral  from  the  fluid  to  the  interface 
and  causing  a  further  growth  of  the  height  of  the  protrusion,  in  a  positive  feedback 
mechcinism.  Selected  wavelengths  grow  exponentially  in  time  [42]. 
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Figure  2-8:  In  reaction  limited  kinetics,  growth  occurs  normal  to  the  interface.  The 
position  of  the  interface  at  time  t  is  indicated  by  a  thick  solid  line  and  the  position 
at  time  t  -h  by  a  dashed  line.  The  change  in  local  height,  Sh{x,t),  relates  to  the 
normal  growth  velocity,  A,  via  Sh^  —  {XSty  +  {XStVhy.  This  relationship  can  be 
derived  from  trigonometrical  considerations,  where  identical  angles  are  indicated  by 
a  double  solid  line  [29]. 
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Figure  2-9;  A  schematic  drawing  of  the  discrete-particle  deposition  “single-step” 
model  of  a  growing  interface  embedded  in  two-dimensional  physiccil  space,  using  pe¬ 
riodic  boundary  conditions.  A  block  of  length  2  is  added  at  each  successive  time  step 
to  a  site  randomly  chosen  among  zdl  sites  that  are  a  local  minima. 
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Figure  2-10:  (a)  Logarithm  (base  10)  of  saturated  value  of  interface  width,  W,  ver¬ 
sus  the  logarithm  of  system  size  L.  For  self- affine  interfaces  W,  will  have  a  linear 
dependence  on  L  with  a  slope  equal  to  3  —  D.  Different  symbols  correspond  to  dif¬ 
ferent  0.5  <  p+  <  1.  Solid  lines  indicate  best  fits  from  linecir  regression.  As 
is  decreased,  the  slope  of  the  graph  decreases,  indicating  an  increase  in  the  fractal 
dimension  of  resulting  interfaces  with  increasing  dissolution- deposition  ratio.  When 
dissolution  begins  to  balance  deposition  {p+  <  0.6),  the  deviations  from  a  linear  de¬ 
pendence  indicate  that  the  interfaces  are  not  self-ciffine.  Error  bars,  taken  as  the 
standard  deviation  of  the  measured  width  from  its  steady  saturation  Veilue,  are  of  the 
order  of  the  symbol  size,  (b)  A  plot  of  W/  versus  ln(L).  Solid  lines  are  best  fits  from 
linear  regression.  This  plot  demonstrates  that  a  better  model  to  fit  the  data  for  high 
dissolution-deposition  ratios  (p+;S0.6),  is  ~  (Ini).  This  transition  is  expected  to 
occur  when  |A|  — »  0  in  2.13. 
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Figure  2-11:  Using  the  slopes  of  the  linear  regression  best-fits  in  Figure2-10a,  we  plot 
D  versus  p./p+j  the  dissolution-deposition  ratio  rate.  Error  bars  are  calculated  using 
the  standard  deviation  from  the  linear  regression  best-fits.  The  fractal  dimension  of 
resulting  interfaces  is  seen  to  increase  with  p_/p+.  For  p_/p+  =  0  D  =  2.63  ±  0.005. 
As  p_/p+  increases  D  approaches  a  limiting  value  of  2.84  ±  0.01.  A  further  increase 
in  p_/p+,  above  p-lp+  ~  0.7,  results  in  suspected  loss  of  fractal  behavior  and  a 
transition  to  logarithmic  dependence,  instead  of  power-law  dependence,  of  W,  on  L. 
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Figure  2-12:  Graphs  of  logio(VF,(L))  as  function  of  logio(-f')  for  (a)  kT  =  \  and  (b) 
kT  =  100.  Solid  lines  are  best  fits  from  linear  regression.  For  kT  =  1  the  fractal 
dimension  increases  (the  slope  decreases)  with  decreasing  similarly  to  results  of 
Model  I.  For  kT  =  100  the  behavior  is  different:  the  slope,  and  hence  the  fractal 
dimension,  is  nearly  constant  with  p_/p+. 


Figure  2-13:  Fractal  dimension  of  interfaces  as  a  function  of  P-jp+,  the  dissolution- 
deposition  ratio.  The  four  curves  correspond  to  four  different  temperatures:  kT  = 
10“^,  10°,  and  10^  from  Model  II  and  kT  =  0  replotted  from  Model  I.  D  is  deduced 
from  the  slopes  of  plots  similar  to  Figures  2-10  and  2-12.  At  high  temperatures  {kT  = 
10^)  the  fractal  dimension  does  not  change  perceptibly  as  a  function  of  p_/p^.  For  low 
temperatures,  curves  follow  a  trend  of  increasing  D  with  increased  p_lp^.  Although 
we  plot  results  from  all  simulations,  interfaces  where  dissolution  and  precipitation 
rates  are  nearly  balanced,  when  p-/p+  >  0.7,  are  probably  not  fractal,  as  can  be 
observed  in  Figures  2-10  and  2-12. 
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Figure  2-14:  A,  the  “roughness  amplitude”,  calculated  from  intercepts  of  lines  in 
Figures  2-lOa  and  2-12a,b  with  the  width  axis,  versus  the  non-dimensioncd  parameter 
P-Ip+.  The  data  follows  a  linear  increase  of  roughness  with  increased  dissolution 
independent  of  kT,  except  for  where  interfaces  are  suspected  not  to  be  fractal.  For 
low  temperatures  when  p.lp+>0.7,  the  roughness  amplitude  A  reaches  a  plateau 
and  diverges  from  the  linear  trend.  Best  fit  to  the  linear  trend  gives  A  =  (0.355  ± 
0.005)p_/p+. 
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Figure  2-15:  For  interfaces  formed  at  low  temperatures  A,  the  roughness  amplitude, 
is  an  increasing  function  of  D,  the  fractcil  dimension,  as  can  be  seen  for  simulation 
results  for  Model  I  (plotted  here  eis  triangles).  Thus,  when  viewed  on  short  length 
scales,  interfaces  with  high  fracted  dimensions  appear  rougher  cind  more  strongly  fluc¬ 
tuating  then  interfaces  with  lower  fractal  dimensions.  Best  fit  to  the  linear  trend  gives 
A  =  (1.01  ±  0.02)D.  The  solid  line  is  theoretical  prediction  given  by  equation  A. 11 
using  values  of  L  =  10^.  Deviation  of  simulation  results  from  theoretical  predictions 
are  probably  due  to  existence  of  natural  cutoffs  for  fractal  behavior,  which  are  not 
represented  in  theory. 
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Measured  Permeability  (mD) 


(B) 


Figure  2-16;  After  [85].  (A)  Demonstrates  the  lack  of  correlation  between  permeabil¬ 
ity  of  the  samples  and  their  measured  total  porosity  squared,  which  is  conven¬ 
tionally  used  for  permeability  predictions  [8].  At  low  porosities  permeability  changes 
more  then  5  orders  of  magnitude  while  is  virtually  constant.  In  contrast,  correla¬ 
tion  with  the  Euclidean  porosity,  is  substantially  better,  since  coincides  with 
the  void  space  relevant  for  transport  calculations.  (B)  Measured  permeability  versus 
permeability  predicted  from  equation  2.18.  The  data  cover  7  orders  of  magnitude  in 
permeability. 
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Chapter  3 


Channeling  instability  of 
upwelling  melt  in  the  mantle 


Abstract 

We  present  results  of  a  theoretical  study  aimed  at  understcinding  melt  extraction  from 
the  upper  mantle.  Specificjilly,  we  address  mechanisms  for  focusing  of  porous  flow  of 
melt  into  conduits  beneath  mid-ocecuu  ridges  in  order  to  explain  the  observation  that 
most  oceanic  residual  peridotites  are  not  in  equilibrium  with  mid-ocean  ridge  basalt. 
The  existence  of  such  conduits  might  also  explain  geological  features,  termed  replacive 
dunites,  that  are  observed  in  exposed  mantle  sections.  We  show  here,  by  linear 
analysis,  that  flow  in  a  chemically  reactive  porous  media  is  unstable  in  the  presence 
of  a  solubility  gradient,  such  as  induced  by  adiabatic  ascent  of  melt  underneath  mid¬ 
ocean  ridges.  The  initially  homogeneous  flow  becomes  focused  in  time  to  produce 
elongated  high-porosity  Angers  that  act  eis  conduits  for  transport  of  fast  flowing  melt. 
This  instability  arises  due  to  a  positive  feedback  mechanism  in  which  a  region  of 
slightly  higher  than  average  porosity  causes  increased  influx  of  unsaturated  flow, 
leading  to  increased  dissolution  which  further  reduces  the  porosity.  Even  in  the 
presence  of  matrix  compaction  emd  chemical  diffusion  the  instability  is  demonstrated 
to  be  robust.  Our  anadysis  also  indicates  the  existence  of  growing,  traveling  waves 
which  transport  and  amplify  porosity  and  concentration  perturbations. 


3.1  Introduction 

Recent  work  [37,  68,  28,  27]  indicates  that  upwelling  mid-ocean  ridges  basalt  (MORE) 
is  in  chemical  disequilibrium  with  the  upper  mantle  peridotites  that  constitute  the 
matrix  through  which  it  flows.  These  observations  place  constraints  on  melting  and 
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melt  extraction  processes  at  ridges.  In  order  to  produce  disequilibrium  transport, 
small  melt  fractions  must  be  efficiently  segregated  from  their  source  and  transported 
to  the  crust  [28,  27,  75,  26].  Since  diffuse  porous  flow  of  melt  along  grain  bound¬ 
aries  would  lead  to  extensive  chemiccd  reaction  and  erasure  of  observed  trace  element 
fractionation,  some  form  of  focused  flow  of  melt  into  channels  has  been  proposed  to 
explain  extraction  of  MORE  from  the  mantle  [79,  22]. 

The  results  of  our  study  imply  that  one  of  the  mechanisms  responsible  for  focusing 
may  be  a  coupled  chemical-hydrodynamical  instability;  uniform  upwelling  of  melt 
flowing  through  a  porous  media  is  unstable  when  the  melt  is  dissolving  some  of 
the  matrix  through  which  it  is  flowing  cind  begins  to  form  elongated,  high-porosity 
channels.  In  ophiolites,  geological  observation  of  dissolution  channels  (dunites)  in 
chemical  equilibrium  with  MORE,  surrounded  by  mantle  peridotites  which  are  not  in 
chemical  equilibrium  with  MORE,  confirms  that  this  instability  may  operate  during 
melt  extraction  from  the  mantle  at  ocecinic  spreading  ridges  [33,  31,  36]. 

Additional  mechanisms  for  melt  extraction  from  the  mantle  beneath  mid-ocean 
ridges  could  include  (1)  hydrofracture  (e.g.,  [56]);  (2)  focused  flow  of  melt  in  zones 
of  localized,  active  deformation  (e.g.,  [81,  35]);  and  (3)  decompaction  into  melt-fllled 
lenses  or  veins  [73].  Mechanisms  1  and  2  eire  most  probable  near  and  above  the  brit¬ 
tle/ductile  transition  in  the  mantle,  where  strain  becomes  localized  into  shear  zones. 
This  is  supported  by  geologiccd  evidence  that  dikes  and  localized  shecir  zones  in  the 
mantle  section  of  ophiolites  form  mostly  “off-axis,”  away  from  a  spreading  ridge,  near 
the  brittle-ductile  treinsition,  and  not  in  the  adiabatically  ascending,  partially  melting 
mantle  beneath  a  spreading  ridge  [35,  32].  The  third  possible  mechanism  is  poorly 
understood  at  present,  and  we  are  not  aware  of  geological  evidence  supporting  such 
a  hypothesis.  In  contrast,  the  reactive  infiltration  instability  is  likely  to  operate  in 
adiabatically  upwelling,  partially  molten,  ductile  asthenosphere,  and  there  is  geologic 
evidence  for  focused  flow  of  melt  in  porous  dissolution  channels  in  the  mantle  section 
of  ophiolites. 

It  has  been  known  for  some  time  [14,  58,  23]  that  reactive  flow  through  a  soluble 
porous  matrix  may  result  in  formation  of  finger-like  embayments  along  an  advancing 
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reaction  front.  The  mechanism  for  this  instability,  termed  the  “reactive  infiltration 
instability”  (RII),  is  simple:  when  unsaturated  fluid  flows  through  a  soluble  matrix,  a 
region  with  slightly  higher  them  average  porosity  will  tend  to  have  an  increased  influx 
of  fluid,  which  will  increase  the  rate  of  dissolution  and  so  increeise  the  porosity  even 
further,  in  a  positive  feedback  mechanism.  Increased  velocity  in  loccilized  regions  will 
cause  lateral  convergence  of  fluid  upstrezun  of  the  front  into  the  high-porosity  fingers 
[36], 

The  characteristic  wavelengths  and  growth  rates  of  the  front  instability  are  deter¬ 
mined  by  three  peirameters:  chemical  reaction  rate,  transport  rate,  and  diffusion  rate 
[80].  Reaction  and  diffusion  act  to  restore  the  system  to  equilibrium,  while  advection 
acts  to  make  it  unstable.  When  diffusion  is  strong  (compared  with  reaction),  it  will 
act  as  the  main  stabilizing  mechanism  for  this  instability.  The  most  unstable  wave¬ 
length  in  this  case  is  determined  by  competition  between  advection  (which  drives  the 
instability)  and  diffusion  (which  tends  to  smooth  perturbations).  The  ratio  between 
these  parameters  is  termed  Pe,  the  Peclet  number. 

Work  on  the  RII  with  reaction-controlled  smoothing  is  sparse  but  generally  indi¬ 
cates  that  growing  fingers  are  present  [24,  80].  Da,  the  Damkohler  number,  is  the 
ratio  between  reaction  timesceile  and  advection  timescale  and  is  the  control  parame¬ 
ter  in  this  case.  Hoefner  and  Fogler  performed  experiments  and  network  simulations 
which  indicate  a  dependence  of  coalescing  or  branching  of  dissolution  channels  on  Da. 

Past  work  on  this  subject  is  not  directly  applicable  to  Earth’s  mantle.  The  front 
problem  as  reviewed  above  is  transient;  there  is  no  supply  of  new  grains  to  the  system, 
and  once  the  reaction  front  has  propagated  through  the  matrix  there  are  no  porosity 
perturbations  left.  Moreover,  the  instability  (area  of  disequilibrium)  is  localized  to 
a  single  interface  between  two  areas  of  equilibrium  rather  than  affecting  the  entire 
interior  solution  (although  [23]  consider  the  case  of  a  front  of  finite  width).  We  see 
no  evidence  for  the  existence  of  such  a  propagating  reaction  front  in  the  mantle,  no 
evidence  of  a  sharp  reaction  zone  underneath  which  the  mantle  is  composed  purely 
of  olivine  eind  above  which  it  is  composed  of  pyroxene.  In  this  paper,  we  investigate 
instabilities  arising  in  a  steady  state  mantle,  where  some  background  porosity,  solid 
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and  liquid  velocities,  and  nndneral  composition  can  be  assumed.  In  this  caise  there  axe 
no  treinsient  solution  fronts,  and  if  an  instability  arises,  it  will  be  present  in  the  whole 
region  of  upwelling  and  dissolution. 

In  the  mantle,  decompression  of  ascending  melt  imdemeath  mid-ocean  ridges 
causes  an  increase  in  solubility  of  solid  phases  with  height  [72,  36].  Choosing  a  point 
along  the  ascent  path  of  the  melt,  one  can  see  that  dissolution  at  this  point  increases 
the  local  concentration  of  soluble  components  in  the  melt  but  never  to  the  point  of 
equilibrium,  since  upward  flow  keeps  bringing  in  undersaturated  melt.  This  small 
departure  from  equilibrium  allows  an  instability  to  occur,  in  the  same  manner  as  the 
feedback  mecheinism  for  the  dissolution  front  described  above.  However,  instead  of 
having  a  Angered  front,  we  expect  the  instability  to  occur  within  the  region  of  melt 
transport  wherever  there  exists  a  gradient  in  solubility. 

In  this  paper  we  study  a  porous  matrix  confined  in  a  box  where  grains  are  sol¬ 
uble,  and  there  is  a  constant  flux  of  melt  from  the  lower  end  of  the  box.  As  solid 
material  dissolves,  the  matrix  is  allowed  to  contract  by  compaction  so  as  to  keep 
the  porosity  constant  in  the  steady  state,  with  additional  grains  supplied  at  the  top 
of  the  box.  The  solubility  of  the  grains  increeises  linecirly  with  increasing  height  in 
order  to  approximate  the  increase  in  solubility  induced  by  adiabatic  ascent  of  melt 
decompressing  in  the  mantle.  Thermcil  melting  of  the  solid  phases,  as  distinguished 
from  reactive  dissolution,  is  neglected,  as  is  viscous  shearing  of  the  solid  phases  and 
advective  heat  transport  by  melt.  In  what  follows,  we  present  the  governing  equa¬ 
tions,  nondimensionaJize  them,  find  a  possible  steady  state,  and  do  a  linear  stability 
eineilysis. 

Two  interesting  unstable  features  are  then  shown  to  coexist: 

1.  The  system  is  shown  to  be  linearly  unstable  to  small  perturbations,  causing 
focusing  of  flow  in  elongated  high-porosity  chainnels,  where  the  vertical  dimension  is 
generally  much  longer  than  the  horizontal  dimension,  establishing  conduits  for  ascent 
of  melt.  These  channels  form  provided  that  the  characteristic  length  for  chemical 
equilibration  is  smaller  than  a  characteristic  length  for  compaction.  A  calculation 
using  characteristic  values  for  Earth’s  mcuatle  predicts  that  the  condition  for  formation 
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of  the  instability  is  probably  met  and  that  the  reaction-infiltration  instability  may 
play  an  importemt  role  in  forming  conduits  for  melt  extraction  from  the  mantle. 

2.  The  system  gives  rise  to  unstable  propagating  waves,  which  in  the  limit  of 
no  dissolution  are  linear  compaction  waves  [76].  The  addition  of  dissolution  during 
porous  flow  gives  rise  to  waves  whose  amplitude  increeises  with  time,  providing  disso¬ 
lution  features  which  propagate  in  space  with  a  finite  phase  and  group  velocity.  These 
results  suggest  a  mechanism  for  spontaneous  nucleation  of  “magmons”  [71,  78]. 

Finally,  we  discuss  the  application  of  our  study  to  focusing  of  melt  flow  in  the 
mantle  beneath  mid-ocean  ridges. 


3.2  Formulation  of  the  Problem 

3.2.1  General  Equations 

In  this  section  we  present  a  set  of  equations  describing  the  essence  of  reactive  flow 
through  a  soluble  porous  medium  with  a  gradient  in  solubility.  The  setup  of  the 
problem  is  given  in  Figure  3-1.  Fluid  is  driven  upward  by  a  pressure  gradient,  entering 
the  soluble  matrix  at  z  =  0  and  leaving  it  at  z  =  L.  With  decreasing  pressure  the  fluid 
has  increasing  ability  to  dissolve  the  porous  matrix.  Since  we  would  like  to  investigate 
a  steady  state  and  deviations  from  it,  we  have  allowed  for  compaction,  though  it  is  by 
no  means  crucial  for  the  growth  of  the  instability.  Compaction  provides  a  relatively 
simple  steady  state  where  dissolution  increases  the  porosity  and  compaction  works  to 
decrease  it. 

The  set  of  governing  equations  presented  below  closely  follows  the  notation  and 
form  of  some  previous  work  on  deformable  porous  media  [44,  71,  76,  77].  This  ap¬ 
proach  views  the  coupled  solid-fluid  system  as  two  interpenetrating  fluids  with  vastly 
different  viscosities  and  is  valid  for  length  scales  much  larger  than  a  pore  size.  Inertial 
effects  have  been  assumed  to  be  negligible. 

Conservation  of  mass.  Conservation  of  the  solid  phase  is  given  by 

^^%^  +  v-i*v(i-^)]  =  -i;ri,  (3.1) 

at  i 
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where  ^  is  the  porosity,  p,  is  the  solid  density  in  kg  m“^,  V  is  the  solid  grain  velocity 
vector  and  Fj  is  the  mciss  trajisfer  rate  of  mineral  i  from  solid  to  fluid  in  kg  m“^  s”l. 
Conservation  of  fluid  meiss  is  given  by 

^  +  v-(/>/v^)  =  Er.-,  (3.2) 

ft 

where  pf  is  the  density  of  the  fluid  and  v  is  the  fluid  velocity  vector. 

Component  conservation  equations  in  the  fluid  phase  consist  of  three  contribu¬ 
tions;  diffusion,  advection,  and  a  chemical  source/sink  term: 

+  V  •  Mvci)  =  V  •  {DiPf4>Vci)  -f  Vi,  (3.3) 

where  D,  is  the  diffusion  coefficient  of  component  i  in  the  fluid  and  c,-  is  the  meiss 
fraction  of  dissolved  component  i  in  the  fluid,  with  2,  c,-  =  1. 

Each  individual  component  is  also  conserved  in  the  solid  phase  such  that 
dp,{l-  <f>)<fi  ^  ^  ^  ^ 

where  cf  and  D'  are  the  mass  fraction  and  the  diffusion  coefficient  of  component  i  in 
the  solid  phase  and  =  1. 

If  one  defines  a  partition  coefficient  Ki  =  c^J Ci  and  assumes  chemical  equilibrium 
between  the  solid  and  the  fluid  phases,  then  previous  formulations  can  be  rederived 
(e.g.,  [44]).  However,  since  we  are  interested  in  nonequilibrium  chemical  reactions, 
we  shall  not  follow  that  practice. 

Mass  transfer  by  chemical  reaction.  Assuming  first-order  chemical  reaction, 
one  can  write  the  rate  of  mass  transfer  as 


Fj  =  -RiA{x,t)  [ci  -  c,,j(j)] , 


(3.5) 


where  Ri  is  the  reaction  rate  constant  of  component  i  in  kg  m~^  s“^ ,  A  is  the  specific 
area  (m*/m^)  available  for  reaction,  cind  Ceq^(z)  is  the  equilibrium  concentration  of 
mineral  i  in  the  fluid  given  in  mass  fraction. 

Solubility  is  teiken  to  be  a  linear  function  of  height,  as  approximately  the  case  for 
melt  that  is  adiabatically  rising  [72,  36]: 


dz 


(3.6) 
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where  Pi  is  a  proportionality  coefficient  describing  the  steepness  of  the  solubility 
gradient  with  units  of 

Darcy’s  law.  Darcy’s  law  relates  the  pressure  p  and  the  relative  velocity  between 
the  fluid  Eind  the  solid  matrix: 


^(v-V)  =  -^Vp,  (3.7) 

fi 

where  the  permeability  of  the  porous  medium  is  usuedly  taken  to  be  a  power  law 
function  of  the  porosity  with  d  a  typical  grain  size,  n  between  2  and  3, 

and  b  a  constant  (e.g.,  [8,  88]).  p  is  the  viscosity  of  the  fluid  and  p  is  the  pressure  in 
excess  of  hydrodynamic  pressure. 

Matrix  deformation.  The  closing  equation  is  the  momentum  conservation  equa¬ 
tion  which  relates  pressure  changes  to  the  rate  of  compaction,  viscous  deformation  of 
the  solid  phcise  and  body  forces  acting  on  the  system  [44,  76]: 


dp 

dxi 


(?Yl  ^ 

dxj  \dxj  dxi ) 


+  •  V  -  (1  - 

oxj  3 


(3.8) 


where  t/,  (  are  the  solid  phase  shear  and  bulk  viscosities  and  Ap  =  />,  -  pf  is  the 
buoyancy  difference.  Equation  3.8  states  that  any  change  in  the  pressure  field  can 
be  expressed  as  the  force  the  solid  exerts  on  the  fluid.  In  a  rigid  material  where  the 
grciin  velocity  goes  to  zero  (V  — >  0),  the  viscosity  of  the  matrix  will  go  to  infinity 
oo)  and  the  product  will  cilways  remain  bounded. 

The  resulting  set  of  dynamical  equations  are  similar  to  those  introduced  by  work¬ 
ers  on  compaction  of  molten  rocks  [44,  76],  but  here  a  specification  of  a  dissolution 
mechcinism  (equation  3.5)  with  increasing  dissolution  as  a  function  of  height  (equa¬ 
tion  3.6)  brings  in  interesting  behavior.  Our  goal  is  to  study  the  combined  effect  of 
dissolution  ^lnd  porous  flow,  rather  than  to  concentrate  on  the  compaction  effects. 


3.2.2  Simplified  Equations 

For  simplicity,  we  assume  the  existence  of  a  fully  soluble  solid  phase  composed  solely 
of  one  mineral  (o'  =  1)  which  can  chemically  react  with  the  fluid  by  dissolution  or 
precipitation,  with  first-order  kinetics.  Since  only  one  reacting  component  is  present. 
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the  subscript  i  will  be  dropped  from  here  on.  The  fluid  phase  is  composed  of  a 
carrier  fluid  with  mass  fraction  1  -  c.  The  carrier  fluid  component  does  not  enter 
the  solid  phase.  The  dissolved  mineral  has  mass  fraction  c  and  effective  reaction  rate 
-Sgff  JiA,  where  we  assume  that  reactive  surface  area  can  be  taken  as  constant  to 
leading  orders. 

The  density  of  the  fluid  phase  is  presumed  constant  as  the  composition  of  the  melt 
changes  due  to  chemical  reaction.  This  assumption  is  approximately  correct  during 
dissolution  reactions  involving  basaltic  melt  and  mantle  minerals.  The  solid  phase 
density  in  a  one  component  system  is  also  a  constant. 

Neglecting  matrix  shear,  the  momentum  conservation  equation  3.8  is 

=  U  +  3^)^^  -  (1  -  4>)^pgk  (3.9) 

where  we  have  defined  a  compaction  rate  as  C  =  V  •  F  and  A  is  a  unit  vector  in  the 
vertical  direction. 

Equations  3. 1-3.5  and  3.7  can  now  be  written  iis 

=  (1  -  4>)C  -  i2gflp(c  -  ceq(z))/p„  (3.10) 

=  --^effCc  -  CeqW)/^/,  (3.11) 

=  DV  .  {<i>Vc)  -  (1  -  c)i2gff(c  -  c^^{z))Ip},  (3.12) 

""  ^  >  (3-13) 

where  3.12  is  a  result  of  subtracting  3.2  from  3.3,  equation  3.4  is  identical  to  3.1  in 
the  case  of  c'  =  1  and  equation  3.13  is  a  result  of  substituting  3.9  in  3.7. 

Boundary  conditions.  In  general,  equations  3.10-3.13  will  require  five  boundary 
conditions  to  solve  for  v,  C,  <f>,  and  c.  Mass  conservation  requires  that  the  flux 
across  a  boundary  be  continuous  or  bcilanced  by  a  source  or  a  sink.  Fliix  boundary 
conditions  include  impermeable,  rigid,  or  “free  flux”  boundary  conditions.  When  the 
boundary  is  impermeable,  the  normal  flux  is  zero  either  because  =  0  or  because 
Vp  •  n  =  0,  where  n  is  the  direction  normal  to  the  boundary.  The  latter  condition 
poses  constraints  on  VC  •  n.  At  a  rigid  boundary  C  =  0,  and  so  in  the  absence  of 


^  +  V.v^ 

f +  V.(V« 
,dc 

<P-^  +  pv  •  Vc 
-^(v  -  V) 
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dissolution  and  matrix  flow,  a  rigid  boundary  has  also  a  constant  porosity  (equation 
3.10).  Finally,  a  free-flux  boundary  has  VC  -n  =  0  and  yields  no  resistance  to  volume 
changes  for  the  normal  flux.  In  addition  to  the  total  flux  boundciry  conditions,  mass 
conservation  poses  constraints  on  flux  of  individueJ  components  in  the  fluid.  These 
constrciints  translate  to  specification  of  concentration  of  solute  in  incoming  or  outgoing 
fluid. 

Nondimensionalization.  In  nondimensioncdizing  equations  3.10-3.13  we  shall 
use  the  following  definitions:  Assuming  zero  solubility  at  z  =  0,  from  3.6, 

ceq(z)  =  /3z. 

The  change  in  solubility  of  the  matrix  between  the  bottom  and  top  of  a  box  of  size 
L  is  defined  as  cs  =  ceq(T)  —  ceq(O),  and  so 

cs  =  L/3.  (3.14) 

The  only  imposed  parameter  value  in  3.14  is  /3,  the  solubility  gradient,  known  from 
thermodynamic  calculations  to  be  of  the  order  of  2  X  10“®  m“^  (an  account  of  calcu¬ 
lations  made  by  [36]  is  given  in  Appendix  B).  If  we  choose  to  investigate  a  small-scale 
box,  then  cg,  the  chcinge  in  solubility  across  the  box,  will  be  small  as  well;  cg  w  0.2 
over  the  vertical  extant  of  decompression  melting,  roughly  the  upper  75  km  of  Earth’s 
mantle  beneath  oceanic  spreading  ridges.  Correspondingly  over  100  m,  cg  is  of  order 
10-^ 

Porosity  is  nondimensionalized  to  a  characteristic  value,  (f>o,  of  the  order  of  10“® 
in  the  partially  molten  upper  mantle  beneath  spreading  ridges  (e.g.,  [27,  75]), 

<f>  =  M'  ■ 


Permeability  is  non-dimensionalized  to  this  porosity 


ko  = 


Fluid  velocity  is  characterized  to  be  of  the  order  of  velocity  driven  by  gravity  forces 
such  that 


^oWo  =  —^P9- 


65 


Chciracteristic  solid  velocity  is  taken  as 


Vo  =  <l>oWoCs. 

Finally,  we  define  a  compaction  length  h  to  be 

k-‘  =  +  4,/3), 

where  >  0  is  an  infinitely  wecik  matrix,  which  compacts  instantly,  and  h  oo  is  a. 
rigid  medium;  h  has  a  typical  value  of  100-1000  m  in  the  mantle  [44]. 

The  nondimensional  variables  will  be  denoted  by  primed  letters: 


{x,z) 

t 

C 

V 

V 


L(x',z') 

WoCs 

Vo^, 

WqV* 

VqV' 

Csc' 


ceq  =  csz 


Where  the  fact  that  time  is  seeded  to  1/cs  is  a  result  of  our  choice  to  scale  time 
to  the  characteristic  time  for  change  in  porosity  due  to  dissolution.  In  the  limit  of 
no  gradient  of  solubility  (eind  thus  no  dissolution),  cs  — >  0,  and  the  characteristic 
timescale  for  cheinge  in  porosity  due  to  dissolution  goes  to  infinity. 

We  now  write  the  nondimensioned  equations: 


^  +  A,V'.V^'  = 

Pt 

(3,15) 

—csDa{c*  —  z'). 

(3.16) 

+  #'v' ■  Vc”  = 

-(1  -  csc')Da{c>  -  z')  ^  V  .  (^i'Vc'), 

(3.17) 

#'(v'-^„V')  = 

(3.18) 
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where  Da  is  defined  as  advection  time  across  a  box  versus  reaction  time: 


Da  =  (3.19) 

•poWoPf 

Alternatively,  if  we  define  the  “reaction  length,” 

(f>OWoPf 

to  be  the  length  scale  over  which  a  perturbation  in  concentration  will  equilibrate  with 
the  solid  matrix  if  it  is  traveling  at  speed  wq,  then  the  Damkohler  number  is  simply 
the  system  size  in  reaction  lengths.  The  reaction  length  and  compaction  length  are 
the  two  inherent  length  scales  in  this  problem. 

The  Peclet  number  is  defined  as  the  advection  rate  versus  diffusion  rate 


Pe  = 


WqL 


(3.21) 


Finally,  we  define  a  rigidity  parameter  a,  which  is  a  combined  measure  of  the  change 
in  solubility  over  one  compaction  length  and  the  size  of  a  system  L  relative  to  a 


compaction  length. 


h?  h 

a  =  —  X  eg  =  —  X  Ch, 


(3.22) 


where  Ch  =  Ph  is  the  change  in  solubility  over  one  compaction  length.  Since  L  caji 
be  as  small  as  we  choose,  a  is  not  necessarily  small  even  if  c/,  is  small.  When  a  0, 
the  matrix  is  easily  compacted,  and  when  a  —*  oo,  the  matrix  is  effectively  rigid. 

We  then  neglect  all  terms  of  order  (since  (f>o  is  of  order  of  10“^).  The  effects 
of  retaining  terms  of  order  when  <  1  have  been  shown  to  be  small  for  many 
problems  [7,  70,  76].  In  addition,  we  temporarily  neglect  diffusion,  Pe  oo,  so  that 
3.15-3.18  become,  dropping  the  primes, 


^  =  C-Da^(,c-z),  (3.23) 

ot  p, 

=  -V  •  {<l>v)  -  cgDa{c  -  z),  (3.24) 

ot 

Cg<f>^  =  -^V  •  Vc  -  (1  -  cgc)Da{c  -  z),  (3.25) 

ot 

=  -<fP[aVC-k].  (3.26) 


67 


The  value  of  cg  increcises  with  system  size,  but  since  we  are  interested  in  physical 
systems  that  have  an  upper  limit  in  size  (the  whole  region  of  decompression  melting 
«  75  km),  Cg  will  be  less  than  or  equal  to  the  concentration  chcinge  over  that  whole 
region,  cg  <  0(10”^),  and  so  will  be  taken  here  as  a  small  enough  parameter  to  allow 
for  expcinsion  techniques. 

Equation  3.23  tells  us  that  the  important  timescale  in  the  problem  is  the  timescale 
over  which  porosity  changes.  This  happens  due  to  compaction  on  the  one  hand  and 
reaction  on  the  other.  Equation  3.24  predicts  that  the  timescale  for  divergence  of 
flux  is  feist  compared  to  that  of  chcinging  the  porosity,  since  cg  is  a  smcill  pcirameter. 
Equation  3.25  predicts  that  the  concentration  in  the  fluid  is  nearly  constant  with  time. 
Finally,  3.26  tells  us  that  pressure  gradients  will  manifest  themselves  as  gradients  in 
compaction  rates. 


3.3  Steady  State 

We  seek  unidirectionail  steady  solutions  to  3.23-3.26,  of  the  form 
[^(z),w{z),c{z),C(z)]  = 

[<t>\z),w\z),c°{z\C\z)]^-cg[<i>\z),w\z),c^{z),C^{z)]^  (3.27) 

where  w{z)  is  the  fluid  velocity  in  the  z  direction,  w°{z)  is  the  zeroth-order  solution 
of  the  steady  state,  eind  cgw^(z)  is  a  small  parameter  correction  to  it;  cg  will  be  shown 
to  be  unimportant  in  the  initicil  steady  state  solution  but  is  included  here  in  order  to 
simplify  the  subsequent  linear  stability  analysis. 

Equations  3.23-3.26  will  be  solved  with  a  physical  picture  in  mind:  at  the  bottom 
of  the  melt  column,  where  melt  is  entering  in  equilibrium  with  its  surroundings, 
there  is  no  dissolution  and  porosity  is  constcint.  Thus  we  require  a  “rigid  boundciry” 
condition  C  =  0,  and  a  chemical  constraint  on  the  concentration  field  at  z  =  0. 
These  lead  to  the  desired  constant  porosity  at  z  =  0.  We  also  require  that  fluid  flux 
is  continuous  across  this  boundary.  Teiking  into  account  that  only  three  boundary 
conditions  are  needed,  now  that  diffusion  of  solute  and  divergence  of  porosity  have 
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been  neglected,  the  boundary  conditions  take  the  form  of: 

^(0)  =  1,  u;(0)  =  l,  c(0)  =  0.  (3.28) 

Solutions  for  the  zeroth-order  flux  are  obtciined  from  the  fluid  conservation  equation 
3.24  and  the  boundary  conditions  3.28 

vPf  =  1,  (3.29) 

From  the  solute  concentration  equation  3.25  one  then  finds  that 

=  +  (3.30) 

where  If  Da  is  a  measure  of  the  thickness  of  a  boundary  layer,  i.e.,  the  dimensionless 
reaction  length  (l/£)a  =  Ireq/f')- 

One  can  see  from  3.30  that  nowhere  (except  for  z  =  0)  does  c°(z)  =  z,  the 
equilibrium  value.  Rather  there  is  always  a  deviation  of  the  concentration  from  the 
equilibrium  concentration,  and  when  the  nondimensional  height  z  »  LeqfL,  the 
deviation  from  equilibrium  approaches  a  constant  undersaturation  (Figure  3-2).  At 
any  point  along  the  ascent  path  of  the  melt,  reaction  tends  to  restore  the  system  to 
equilibrium,  but  more  undersaturated  fluid  is  brought  from  below  to  drive  the  system 
away  from  equilibrium.  This  solution  suggests  that  some  degree  of  disequilibrium  will 
exist  eis  long  as  the  reaction  length  is  significantly  larger  than  the  continuum  length 
scale  (i.e.,  the  grain  scale).  When  Leq  <  0{d),  then  the  system  is  effectively  in  local 
equilibrium.  Note  that  we  will  show  that  the  channeling  instability  arises  even  under 
conditions  of  effectively  local  equilibrium. 

For  simplicity,  we  shall  assume  that  the  nondimensional  height  is  z  LeqlL 
such  that  any  boundary  layer  effects  are  negligible.  Estimates  of  the  reaction  length 
that  are  reasonable  for  Earth  are  sensitive  to  assumptions  about  the  microscopic 
distribution  of  melt  and  solid  (see  Appendix  B).  For  a  reinge  of  parameters.  Table  3.1 
suggests  that  the  equilibration  length  may  range  from  much  less  them  a  millimeter 
to  meters.  So  even  if  our  system  size  is  of  the  order  of  a  compaction  length,  most 
parameter  rzmges  indicate  that  significant  boundeiry  layers  are  not  expected.  Thus 
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we  can  approximate 


c°(z)  «  2  (3.31) 

The  zeroth-order  steady  state  compaction  rate  is  calculated  from  3.23  to  be  nearly 
constant  for  x  L&iilL 

C°(z)  =  -  1)  «  (3.32) 

P*  P$ 

From  3.26  the  zero-order  steady  state  porosity  is  approximately  constant: 

e“^“^  4-1^  ss  1.  (3.33) 

Finally,  from  3.29  the  velocity  is  also  approximately  constant: 

w%z)  «  1.  (3.34) 

The  zeroth-order  steady  state  of  constant  porosity  and  upwelHng  velocity  is  sustained 
by  the  competition  between  porosity  formed  by  dissolution  and  destroyed  by  com¬ 
paction.  Since  compaction  rate  is  the  gradient  of  grain  velocity,  constant  compaction 
means  that  grains  descend  with  increasing  speed  as  a  function  of  height,  which  exactly 
balances  the  net  incrccise  in  dissolution  with  height. 

The  C5-order  terms  of  the  steady  state  can  be  easily  obtained  from  3.23-3.26  using 
the  zeroth-order  solutions  and  boundary  conditions  of  ^^(0)  =  lu^(O)  =  c^(0)  =  0,  but 
their  detail  is  of  no  particulcir  interest  here,  since  in  the  stability  analysis  perturbations 
to  terms  of  order  cs  axe  negligibly  small  compared  to  perturbations  to  terms  of  order 
0. 


3.4  Linear  Analysis  without  diffusion 

We  shall  perform  a  linear  stability  analysis  of  3.23-3.26  assuming  that  all  variables 
Ccin  be  expressed  cis  their  steady  state  Veilue  plus  small  deviations: 

t),  w\x,  t),  u\x,  t),  c'(x,  i),  C'(x,  t)]  = 

[^(z),  w(z),  0,  c(z),  C{z)]  -h  e[^(x,  t),  iw(x,  t),  u{x,  t),  c(x,  t),  C(x,  i)],  (3.35) 
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where  e  -C  1  and  the  steady  state  values  are  defined  by  3.27-3.34. 

Keeping  0(e)  terms  and  discarding  terms  of  0{cse,  e^),  the  perturbation  equations 
take  the  following  form: 


=  C-Da^c, 

(3.36) 

Pi 

dw  du 

dz  ^  dz  ^  dx^ 

(3.37) 

A  dc 

=  — ^  —  w  —  — - Dac, 

oz 

(3.38) 

1  lU  ^ 

=  (n-l)#-a^. 

(3.39) 

dC 

“ax’ 

(3.40) 

In  solving  the  perturbation  equations  3.36-3.40  we  shall  assume  that  all  variables 
have  the  form  of 

[u;(x,  i),  u(x,  t),  ^(x,  i),  c(x,  t),  C(x,  t)]  =  [w{z),  u{z),  ^(z),  0(2),  C(z)]e‘"‘e‘'*,  (3.41) 

where  a  is  the  nondimensional  growth  rate  of  the  perturbation  and  I  is  the  nondi- 
mensional  wavenumber  in  the  horizontal  direction.  Equation  3.38  can  be  rewritten 
using  3.36,  3.39,  and  3.40  to  be  a  function  of  C  only 

nDa^  -  a(V  +  Da)  ^  =  {aDa^  -  1)V  -  Da  C,  (3.42) 

Pi  J  L 

where  djdz  is  designated  by  the  operator  symbol  V. 

Equation  3.37  can  also  be  rewritten  using  C, 

nV^  =  a{V^  -  l^)C,  (3.43) 

eliminating  ^  from  both  of  these  equations,  one  arrives  at  a  final  equation  for  a  single 
variable 

+  {aaDa  —  n)V^  —  {l^acr  -1-  nDa)V  -1-  al^Da{npf  f  p,  —  o-)j  C  =  0.  (3.44) 

3.4.1  Preview  of  Solutions  and  a  Simple  Scaling  Argument 

There  cire  two  kinds  of  instabilities  that  we  find  in  this  system:  One  is  the  growth 
of  an  absolute  instability,  which  is  stationary  in  space  and  obeys  a  set  of  boundary 
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conditions.  In  this  Ccise  the  growth  rate  c  is  purely  real  and  so  dissolution  features 
will  grow  pinned  in  space  and  will  not  travel.  The  other  kind  of  instability  is  growth 
in  time  of  the  amplitude  of  traveling  waves  (in  this  case,  a  is  complex),  in  which 
we  assume  a  semi-infinite  medium  and  investigate  the  behavior  of  traveling  waves 
without  imposing  boundary  or  initial  conditions  (these  can  be  imposed  in  a  future 
investigation  in  order  to  study  the  effects  of  the  finite  size  of  the  medium  on  the 
imstable  wave  solutions). 

How  do  we  expect  the  most  unstable  wavelength  of  the  absolute  instability  to 
behave  as  we  change  Da,  the  control  parameter,  and  should  we  expect  a  dominant 
wavelength  to  emerge  at  till?  As  a  perturbation  of  horizontal  wavelength  A*  grows, 
unsaturated  fluid  converges  laterally  toward  growing  features:  when  there  is  no  dif¬ 
fusion  only  reaction  is  present  to  counteract  and  the  deviations  from  equilibrium, 
we  expect  that  if  the  time  to  advect  laterally  across  a  perturbation  (A^/u)  is  longer 
than  time  for  reaction  to  wipe  out  the  concentration  difference  (l/i2g^),  then  the 
perturbation  will  not  focus  enough  unsaturated  fluid  to  keep  itself  alive  and  it  will 
be  damped.  This  means  that  perturbations  with  long  horizontal  wavelengths  com¬ 
pared  to  Xeq  will  not  grow  effectively.  On  the  other  hand,  focusing  by  the  longest  of 
the  fast  growing  wavelengths  will  starve  the  shorter  ones,  and  a  dominant  horizontal 
wavelength  is  expected  to  emerge.  By  this  argument,  the  horizontal  wavelength  of  the 
most  unstable  mode  should  increase  with  increasing  Leq  (decreasing  Da  number). 

Compaction  is  expected  to  damp  horizontal  wavelengths  comparable  to  a  com¬ 
paction  length.  We  propose  that  if  Leq  is  so  large  that  that  the  most  unstable 
wavelength  is  of  the  order  of  a  compaction  length,  then  stationary  channels  could 
not  be  maintained  in  the  system.  However,  the  results  of  our  study  indicate  that 
even  when  stationary  channels  are  inhibited  from  growth,  the  system  still  exhibits  a 
“traveling  instability”:  unstable  growth  of  traveling  waves. 

3.4.2  Unstable  Stationary  Channels 

In  this  section  we  investigate  the  growth  of  unstable  dissolution  features  by  inves- 
tigating  the  growth  rate  of  vertical  modes  that  obey  boundary  conditions,  termed 
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“absolute  instability.” 

A  solution  to  3.44  of  the  form 

C{z)  =  +  Aae"**'  (3.45) 

will  exist,  provided  that 

+  {Da  — —)m^  —  {l^  +  ^^--)rn  +  —  1)  =  0,  (3.46) 

aa  aa  ptO 

where  mi,2,3>  the  three  roots  of  the  cubic  polynomial  3.46,  are  either  all  real  or  one 
real  and  two  complex  conjugates.  Equation  3.46  establishes  the  relationship  be¬ 
tween  wavenumbers  in  the  vertical  (mi, m2, m3)  and  horizontal  {1)  directions  and 
their  growth  rate.  In  order  to  find  the  solution  of  the  linear  stability  problem  (that 
is,  to  find  the  growth  rate  <7  as  a  function  of  the  Damkholer  number  and  the  rigidity 
a  for  any  given  horizontal  wavenumber  /)  we  need  to  specify  a  set  of  three  boundary 
conditions,  which  actually  correspond  to  the  third-order  differential  equation  3.44. 
These  boundary  conditions  will  constrain  the  vertical  modes  mi,  m2  and  m3,  and  as 
a  consequence  determine  a{l,  Da,  a)  via  3.46. 

Boundary  conditions  for  equation  3.44  emerge  from  the  following  assumptions: 
(1)  There  is  no  dissolution  and  porosity  is  constant  at  z  =  0,  where  the  incoming 
fluid  is  in  equilibrium  with  the  matrix.  This  assumption  leads,  via  3.36,  to  a  rigid 
boundciry  condition  (C  =  0)  at  z  =  0.  (2)  The  z  =  0  boundary  is  impermeable 
to  the  perturbation,  meming  that  the  flux  of  fluid  normal  to  the  boundary  remciins 
unperturbed  from  its  steady  state  value.  (3)  Using  an  observation  from  physical  and 
computer  experiments  [36]  that  lateral  fluxes  ahead  of  the  perturbation  are  negligible, 
we  require  (from  3.40)  a  rigid  boundary  at  z  =  Z  as  well.  (Alternatively,  one  could 
choose  a  “free-flux”  boundary  at  z  =  Z,  which  actually  acts  to  amplify  the  pertur¬ 
bation  by  relaxing  the  restrictive  top  rigid-boundary  condition,  and  also  complicates 
the  mathematical  presentation  somewhat.)  The  above  conditions  are  equivalent  to 

dC 

C{z  =  0)  =  0,  C{z  =  1)  =  0,  -^{z  =  0)  =  0,  (3.47) 

and  we  seek  the  conditions  under  which  a  nontrivial  solution  of  the  form  3.45  exists. 
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The  boundciry  conditions  3.47  tell  us  that 


''ll  1  ^ 

(  A  \ 

Ax 

fo) 

gmi  ginj  gmj 

A2 

= 

0 

^  mi  m2  m3  ^ 

This  has  a  nontrivial  solution  if  the  determinant  is  equal  to  zero: 


(3.48) 


-m3)  +  e”‘’(m3-mi)  +  e”**(mi -m2)  =  0.  (3.49) 

To  find  the  growth  rate  as  a  function  of  horizontal  wavenumber  from  3.46  and  3.49,  we 
first  analytically  calculate  the  three  roots  mi^2,3(o‘>  ot)  of  the  cubic  polynomial 
3.46.  Substituting  these  roots  into  3.49,  we  then  obtain  an  implicit  equation  for 
a{l,  Da,  a).  Choosing  a  Vcdue  for  the  parameters  Da  and  a,  we  finally  obtain  £7(1)  by 
seeking  the  roots  of  the  implicit  equation  3.49,  using  a  bisection  numerical  method. 
Numerical  solutions  indicate  that  mi, m2, m3  have  the  form  of  1  real  root  and  2 
complex  conjugates  so  that  C{z)  of  3.45  can  be  written  as 

C{z)  =  5ie"“*  +  cos  qz  +  sin  qz.  (3.50) 


Before  presenting  the  results  of  the  linear  stability  analysis,  we  would  like  to  briefly 
discuss  the  physics  of  the  problem  revealed  by  writing  3.46  as  a  dispersion  relation  in 
which  £7  =  £7(m,  /): 


£7  = 


1  nm 


nPDa^ 

Pi 


am?  —  B  (m*  — /2)(m  +  7?a)  (3.51) 

The  growth  rate  a  is  composed  of  two  completely  separate  parts,  one  that  includes 
a  dependence  on  the  rigidity  of  the  matrix,  a,  but  does  not  depend  on  the  chemical 
reaction,  and  the  second  that  depends  only  on  the  rate  of  chemical  reactions,  the  Da 
number,  but  does  not  depend  on  the  rigidity;  a  can  be  expressed  as  a  sum  of  these 
parts 

a  =(Tc+  (TDa,  (3.52) 

where  the  compaction  contribution  to  the  growth  rate  is 

1  nm 


£7c  = 


a  {rn?  —  P) 


(3.53) 
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and  the  chemical  reaction  contribution  to  the  growth  rate  is 


<^Da  = 


nPDa^ 

_ £» _ 

(m?  —  P){m  +  Da) 


(3.54) 


As  the  matrix  rigidity  is  increased  (a  — >  oo),  then  a  — >  ctdo-  This  limit  strips  away 
compaction  effects  on  the  instability.  In  the  limit  of  Da  —*  0,  a  —*  ac  and  only 
compaction  effects  cire  left. 

Rigid  medium  limit.  As  mentioned  previously,  in  the  limit  of  a  rigid  medium 
(1/a  0),  compaction  effects  are  not  present,  illuminating  the  physics  of  the  dis¬ 

solution  instability.  In  Figure  3-3  we  present  the  growth  rate  cr  as  a  function  of  the 
horizontal  wavenumber  I  for  Da  =  10  in  a  rigid  medium  (1/a  =  0).  For  compari¬ 
son,  the  growth  rate  in  a  compacting  medium  with  a  =  1  is  also  shown.  The  choice 
of  Da  —  10  is  given  as  an  example,  to  demonstrate  the  qualitative  behavior  of  the 
solution.  For  a  rigid  medium  any  value  of  Da  produces  the  same  kind  of  behavior 
with  positive  growth  rate  peaking  at  a  certain  wavelength.  Results  for  a  compacting 
medium  will  be  discussed  in  the  next  section. 

It  should  be  noted  that  actually  cr  attains  several  values  for  each  value  of  1:  these 
correspond  to  different  growth  rates  of  different  vertical  modes  with  a  wavelength 
A,  =  27r/g  where  q  is  defined  in  3.50.  The  first  mode  hcis  approximately  half  a 
wavelength  (q  close  to  tt)  in  the  vertical  dimension  of  the  box  and  is  the  fastest 
growing  mode.  The  second  mode  has  close  to  one  wavelength  fitted  in  the  vertical 
dimension  and  grows  more  slowly.  The  third  mode  grows  even  more  slowly,  etc.  Hence 
only  the  first  mode,  the  fastest  growing  one,  is  plotted  on  Figure  3-3.  The  growth  rate 
in  Figure  3-3  is  seen  to  peak  for  horizontal  wavenumber  Imax,  and  so  Ax  =  27r/Zrnax 
is  the  most  unstable  wavelength  in  the  system. 

Calculations  similar  to  Figure  3-3  have  been  made  for  different  Da  numbers.  Fig¬ 
ure  3-4a  shows  the  most  unstable  wavelengths.  A*,  as  a  function  of  Da,  for  the  rigid 
medium  case.  As  Da  increcises.  A*  is  shown  to  decrecise  and  to  approach  a  power  law 
dependence  on  Da.  For  Da  ^  1,  and  spanning  4  orders  of  magnitude,  the  nondimen- 
sional  dominant  horizontal  wavelength  scales  as  AJ,.  ~  Jl/Da.  In  dimensional  units 
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this  means  that 

Ax  ~  *  ieq.  (3.55) 

The  increase  in  dominant  horizontal  wavelengths  of  chemnels  with  increasing  equili¬ 
bration  length  (decreasing  Da)  is  cis  predicted  by  the  “preview  of  solutions.”  The 
fact  that  the  lateral  extent  of  channels  depends  also  on  the  verticcil  dimension  of  the 
box  is  more  surprising.  We  postulate  that  the  verticeil  length  scale  L  is  imposed  on 
the  perturbations  by  the  fact  that  channels  edways  span  the  box  verticeiUy.  Since  the 
vertical  mode  is  coupled  to  the  horizontal  modes  in  3.46,  the  horizontcil  modes  are 
forced  to  feel  the  system  size  too.  The  aspect  ratio  (Ax/A^)  of  the  channels  decreases 
as  Da  is  increased.  Since  in  dimension<il  parameters  (when  Da  1)  L  and 

Ax  ~  \JL  *  Leq,  then  the  aspect  ratio  is 


as  seen  in  Figure  3-4b,  cind  the  chcinnels  become  more  finger-like  as  the  equilibration 
length  decreases  or  as  the  system  size  increases. 

The  growth  rate  cr  of  the  fastest  growing  horizontal  wavenumber  /max  as  a  function 
of  Da  is  plotted  in  Figure  3-4c,  where  one  can  see  that  as  Da  — >  oo,  the  growth  rate 
approaches  a  constant  Vcdue,  a  — >  np//p„  which  can  be  predicted  from  3.54.  This 
limit  is  determined  by  dkfd<f><x  n,  the  derivative  of  the  permeability  with  respect  to 
porosity.  If  the  permeability  decreased  with  increasing  porosity,  then  the  instability 
would  not  occur  and  the  growth  rate  would  approach  a  negative  constant  value.  In 
other  words,  the  maximum  change  in  porosity  is  only  related  to  the  rate  at  which 
flux  changes  with  porosity. 

The  increase  of  cr  with  Da  seems  counterintuitive  at  first,  since  as  Da  is  increased 
the  system  is  closer  to  equilibrium  (equation  3.31  and  Figure  3-2).  The  explcination 
stems  from  the  fact  that  as  Da  — >  oo,  emy  perturbation  in  flux  is  immediately  com¬ 
pensated  by  chemical  reaction  bringing  the  liquid  close  to  local  chemical  equilibrium. 
In  the  mecintime,  the  porosity  has  been  lowered  further  by  the  strong  dissolution,  so 
additional  fluid  flows  into  the  perturbed  region.  This  in  turn  will  result  in  immediate 
strong  dissolution  lowering  the  porosity  even  further.  On  the  other  hand,  a  low  Da 
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number  will  tend  to  weaken  the  growth  of  channels  because  the  fluid  has  time  to 
redistribute  the  concentration  by  lateral  advection  and  thus  smooth  away  gradients 
before  substcintial  dissolution  occurs. 

Compacting  medium.  Numerical  solutions  for  a  compacting  medium  with  a  of 
order  1  indicate  that  there  exists  a  critical  value  for  Da.  For 

Da  >  Da^rit  =  (3-57) 

dissolution  channels  will  grow.  Once  condition  3.57  is  violated,  there  are  no  solutions 
which  obey  the  boundary  conditions  3.47  and  dissolution  channels  do  not  form. 

Condition  3.57  can  be  rewritten  using  the  definitions  of  Da,  a  and  L  (equations 
3.19,  3.22,  and  3.14): 

Daa  =  >  1.  (3.58) 

Leq 

In  other  words,  this  term  is  proportional  to  the  ratio  of  the  compaction  length  to  the 
reaction  length.  If  the  compaction  length  is  much  larger  than  the  reaction  length, 
permanent  channels  can  grow.  Otherwise,  compaction  will  become  important  and 
will  act  as  a  stabilizing  mechanism,  trying  to  squeeze  the  channels  shut  and  forcing 
propagation  of  the  instability  as  waves,  cis  can  be  seen  in  the  “unstable  dissolution 
waves”  section. 

Figure  3-3  shows  the  growth  rate  for  Da  =  10,  a  =  1,  so  that  Da  a  =  10,  not  so 
high  above  the  critical  value  Jind  in  the  lower  range  of  values  of  Da  X  a  predicted  for  the 
mcintle,  as  can  be  seen  in  Table  3.1.  Growth  rates  in  a  compacting  medium  are  shown 
in  comparison  with  results  for  Da  =  10  and  a  rigid  medium,  1/a  =  0.  Figure  3-3 
shows  that  the  maximum  wavelength  in  the  compacting  medium  is  not  changed  from 
a  completely  rigid  medium  eind  the  growth  rate  is  only  slightly  lower  than  in  the  rigid 
case.  Compaction  is  shown  to  damp  the  long  wavelengths,  as  expected. 

Eigenfunctions.  Here  we  calculate  the  eigenfunctions  (full  z  dependent  solu¬ 
tions)  of  all  perturbed  fields  in  the  problem  {C,^,u,v,c)  and  plot  them  as  contour 
plots,  taking  a  snapshot  in  time.  This  demonstrates  our  claim  that  the  instability  is 
characterized  by  channel-like  features  and  helps  one  to  visueilize  the  spatial  distribu¬ 
tion  of  melt  and  porosity. 
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We  use  3.50,  the  three  independent  boundary  conditions  3.47,  and  3.36-3.41  to 
analytically  find  the  eigenfunctions  and  thus  the  full  solutions  for  the  perturbation 
variables: 


C(x,z,t)  — 


u(x,z,t)  = 
w(x,z,i)  = 

c(x,z,i)  = 

In  order  to  plot  the  behavior  of  the  perturbation  variables  as  a  function  of  x  and 
z,  we  use  the  values  of  mi^p,q,  and  a  calculated  numerically  for  the  most  unstable 
wavelength,  as  described  in  the  previous  section.  In  the  surface  plots  illustrated  in 
Figure  3-5  we  used  parameter  values  of  Da  =  100,  a  =  1. 

Figure  3-5a  is  a  plot  of  the  perturbation  in  compaction  C{x,  z,  to),  which  is  equiv¬ 
alent  to  the  perturbation  pressure.  Narrow  channels  spanning  the  vertical  dimension 
of  the  box  can  be  seen.  The  constant  pressure  boundary  conditions  force  the  pressure 
field  to  achieve  a  maximum  at  about  three  quarters  of  the  way  to  the  top  of  the  box 
and  not  at  the  top  boundary.  Porosity  ^{x,  z,  to)  is  plotted  in  Figure  3-5b  and  can  be 
seen  to  have  increasing  amplitude  with  increasing  height  and  to  achieve  a  maximum 
amplitude  at  the  top  of  the  box.  Disequilibrium  undersaturation,  —c[x,z,to),  and 
vertical  velocity  have  surface  plots  very  similar  in  shape  to  the  plot  of  porosity. 
Within  linear  approximations,  when  the  system  is  unstable,  the  longer  the  melt 
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ascends,  the  more  robust  the  instability  becomes,  spanning  any  given  box  size.  This 
is  understandable,  since  the  more  disequilibrated  the  ascending  fiuid,  the  stronger  the 
driving  force  for  instability. 
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3.4.3  Unstable  Dissolution  Waves 

In  this  section  we  present  a  different  solution  to  the  set  of  perturbation  equations 

3.36- 3.40,  which  is  independent  of  the  channel  solution  and  may  coexist  with  it.  We 
study  simple  linear  propagating  waves  in  an  infinite  medium.  These  are  related  to 
compaction  waves  but  present  a  previously  unstudied  facet  of  the  RII  instability, 
unstable  traveling  dissolution  waves. 

When  mciss  transfer  between  the  solid  and  the  liquid  is  zero  (Da  =  0),  equations 

3.36- 3.40  are  close  to  the  ones  arrived  by  [76]  in  his  linear  analysis  of  compaction 
driven  waves,  apart  from  another  time  derivative  which  survived  in  his  divergence 
equation  because  of  the  different  time  scciles  and  therefore  different  linearization  in¬ 
volved.  In  the  limit  of  £>a  — *  0  we  arrive  at  an  equation  for  traveling  waves: 

n|-V'|  =  0.  (3.64) 

Assuming  wave  solutions  of  the  form 

4>  =  (3.65) 


the  dispersion  relation  is 


nm 


cr  = 


indicating  the  existence  of  traveling  waves,  with  a  phase  velocity  Cp 

cr  n  cos  9 


(3.66) 


K 


K2  ’ 


(3.67) 


where  the  wave  vector  is  defined  by  and  9  is  the  angle  between  the  wave 

vector  and  the  vertical,  as  seen  in  Figure  3-6a.  The  phase  velocity  in  the  z  direction, 
Cp,  =  Cj,f  cos  9,  is  greater  than  zero,  thus  linear  compaction  waves  always  travel  with 
an  upward  component.  (It  is  interesting  to  note  that  3.64  is  identical  to  the  equation 
for  planetary  Rossby  waves  that  dominate  the  ocean  and  atmosphere.  Rossby  waves 
have  westward  traveling  phases.) 

These  results  from  the  zero  Da  limit  are  close  to  the  results  from  [76],  in  which 
traveling  dispersive  waves  arise  from  compaction.  The  waves  form  due  to  the  increase 
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of  melt  flux  as  a  function  of  porosity,  and  its  ability  to  deform  the  matrix.  Viscous 
resistance  to  volume  chcinge  causes  waves  to  disperse. 

We  now  turn  to  study  the  full  problem  of  nonzero  Da  and  seek  wave-like  solutions 
to  the  set  of  perturbation  equations  3.36-3.40, 

C{z,x,t)  ~ 


where  m  and  I  are  reed  and  a  is  to  be  determined  from  the  dispersion  relation.  A 
single  equation  similar  to  3.46  is  found,  which  in  turn  Ccin  be  written  as  a  dispersion 
relation  that  restricts  the  growth  rate  to  be 

m  h  ^ 

(m^  +  D)(m^  -f  Da^)  +  D  a  p,{m^  -f  Da^) 

This  growth  rate  has  an  imaginary  part,  aj,  and  a  reed  positive  part,  a-z,  indicating 
that  if  porosity  waves  were  formed  in  the  presence  of  dissolution,  their  amplitude 
would  increase  with  time  to  form  “unstable  traveling  dissolution  waves”,  or  “chan- 
nelleons” . 

The  phase  velocity  of  the  waves,  Cp,  is 


(3.68) 


/  1  sin^  \ 

\aK^  ^  K^cos2e  + Day  ' 


(3.69) 


where  the  wave  vector  is  defined  by  p  +  and  9  is  the  angle  between  the  wave 
vector  and  the  vertical.  In  the  presence  of  dissolution,  the  amplitude  of  the  waves 
grows  in  time  due  to  the  positive  reed  part  of  the  growth  rate. 


nsin^  9Da^^ 

_ pt 

cos^  6  -f  Da^ 


(3.70) 


Here  an  is  seen  to  be  edways  positive,  providing  a  mechanism  for  nucleation  and 
growth  of  “magmons” ,  which  previously  required  an  initial  step  perturbation  in  poros¬ 
ity  in  order  to  nucleate.  For  Da  =  0,  erz  =  0,  indicating  that  compaction  waves  in 
the  absence  of  dissolution  axe  marginally  stable  and  travel  with  a  constant  amplitude. 

The  presence  of  dissolution  (Da  >  0)  brings  about  interesting  behavior  of  both 
the  phase  velocity  of  the  waves  and  the  rate  at  which  their  amplitude  grows.  Figure  3- 
6b  illustrates  Cp  and  oz  as  calculated  from  3.69  and  3.70.  The  top  panel  shows  Cp 
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and  0-71  as  a  function  of  orientation  of  the  wave  front  9  for  consteint  K,  Da,  and  a. 
The  phase  velocity  of  planar  waves  is  meiximal  when  they  have  a  diagonal  orientation 
but  drops  to  zero  for  waves  that  are  oriented  vertically  {6  — >  Tr/2).  On  the  other 
hand,  these  vertical  “stationary  channels”  are  the  ones  that  grow  the  fastest,  as  is 
seen  from  (t{0),  and  are  actually  the  collapse  of  the  traveling  wave  solution  to  the 
stationary  chcinnel  solution  obtained  in  the  previous  section.  The  bottom  panel  of 
Figure  3-6B  shows  Cp  and  as  a  function  of  Da  for  constant  K,  9,  and  a.  The 
maximum  in  phase  velocity  corresponds  to  K  cos  9  =  Da,  while  maximal  growth 
rate  occurs  for  Da  00,  similar  to  the  maximal  growth  rate  for  the  stationary 
channels.  Compaction-dissolution  waves  are  dispersive,  as  expected  for  waves  which 
arise  partly  from  compaction  [76],  and  the  long  wavelengths  both  travel  and  grow  the 
fastest.  Like  the  linear  compaction  waves  of  3.64,  these  waves  always  have  an  upward 
phase  velocity.  Waves  propagate  with  a  finite  phase  velocity  even  in  the  absence  of 
compaction.  They  appear  to  travel  diagonally  because  vertical  wavelengths  that  are 
shorter  than  the  system  size  (m  >  0  or  5  <  7r/2)  imply  that  at  some  point  in  space, 
upwelling  fluids  in  high-porosity  regions  encounter  an  obstacle  of  low  porosity,  forcing 
the  fluid  to  chew  its  way  up  with  a  diagonal  component  of  velocity. 


3.5  Finite  Diffusion 

Diffusion  in  a  porous  medium  is  a  result  of  hydrodynamic  dispersion  and  molecular 
diffusion.  Dispersion,  in  most  cases,  is  associated  with  larger  diffusion  coefficients  than 
moleculcir  diffusion.  Both  have  the  effect  of  of  causing  an  initial  sharp  concentration 
gradient  to  spread  out  with  time.  This  section  addresses  the  question  of  whether  a 
finite  diffusion  rate  will  smooth  out  gradients  in  the  concentration  field  to  a  point 
where  the  channeling  instability  will  not  be  able  to  grow. 

3.5.1  Predictions 

It  is  probable  that  perturbations  smaller  than  a  diffusion  length  will  be  smoothed 
out.  When  diffusion  is  weak,  and  a  diffusion  length  is  smaller  than  the  most  unstable 
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wavelength  determined  by  reaction  (as  calculated  for  the  no  diffusion  ceise  above), 
diffusion  will  only  act  as  a  short-wave  cutoff.  When  diffusion  is  strong,  and  a  diffusion 
length  is  longer  than  the  most  unstable  wavelength  determined  by  reaction  alone,  it 
will  allow  for  unstable  growth  but  modify  the  dominant  horizontal  wavelength,  forcing 
it  to  become  longer  than  a  diffusion  length.  This  is  because,  as  seen  in  Figure  3-3, 
all  wavelengths  are  unstable  in  the  presence  of  reaction  alone  (in  the  rigid  case)  and 
when  short  waves  are  damped,  long  wavelengths  can  become  the  fastest  growing  in  the 
system.  It  follows  that  for  strong  diffusion,  the  growth  rate  will  be  lowered  compared 
to  the  case  with  no  diffusion. 


3.5.2  Simplified  Calculations 

We  present  here  a  somewhat  simplified  linear  analysis  for  the  finite  diffusion  case.  We 
allow  for  diffusion  only  in  the  horizontal  direction,  since  allowing  for  vertical  diffusion 
adds  higher  order  derivatives  in  the  vertical  direction,  complicating  the  vertical  struc¬ 
ture  of  the  equations  beyond  the  point  necessary  to  obtain  well-behaved  solutions  to 
the  stability  problem.  This  approximation  is  probably  justified,  since  the  vertical 
melt  velocity  is  believed  to  be  much  larger  than  diffusion  rates.  However,  one  should 
note  that  available  experimental  diffusivities  in  silicate  melts,  quoted  in  Table  3.1,  do 
not  include  effects  of  dispersion. 

Solving  the  perturbation  equations  3.36-3.40  with  a  finite  value  for  1/Pe,  3.46  will 
have  the  form 

m"  -f  {Da  +  /VPe  - —)m^  -  (P  -h  — )m  -1-  -  l\Da  -f  P/Pe)  =  0.  (3.71) 

aa  aa  p,cr 

The  growth  rate  is  again  decomposed  into  compaction  and  reaction  parts 


<^  =  Crc  +  O^Da 

(3.72) 

where  the  compaction  contribution  to  the  growth  rate  is 

1  nm 

(3.73) 

^  a  {m}  —  P) 

and  now 

(m2  -  P)(m  +  Da  ^  PjPe) 

(3.74) 
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is  the  diifusion-reaction  contribution  to  the  growth  rate.  In  the  limit  of  /  ->  oo, 
»  0,  as  expected  since  diffusion  will  tend  to  damp  the  growth  of  short  wavelength 
perturbations. 

We  use  boundary  conditions  3.47  to  solve  3.71  following  the  same  numerical  pro¬ 
cedure  outlined  in  the  “unstable  stationary  chcinnels”  section.  Figure  3-7  demon¬ 
strates  the  effect  of  adding  diffusion  to  the  no-diffusion  solution  obtained  for  the  rigid 
medium.  Generally,  diffusion  acts  to  reduce  the  growth  rate  of  the  instability  and  to 
shift  the  most  unstable  wavelength  to  become  longer  than  in  the  no  diffusion  case,  as 
predicted. 

Figure  3-8  shows  the  most  unstable  wavelength  and  its  growth  rate  for  a  constant 
Da  =  10  and  a  varying  Pe  number.  The  qualitative  behavior  is  as  predicted:  for 
weak  diffusion  (Pe  :§>  Pa,  which  corresponds  to  the  condition  when  equation  3.17  can 
be  considered  diffusionless)  the  most  unstable  wavelength  and  its  growth  rate  remain 
nearly  the  same  as  the  case  with  no  diffusion.  For  Pe  <  Da  the  dominant  wavelength 
A,  increases  with  increasing  diffusion  length  (decreasing  Pe)  A,  oc  (l/Pe)^/'*,  and  the 
growth  rate  is  reduced  but  remains  positive. 

Thus  we  conclude  that  when  diffusion  is  strong  or  when  the  Pe  number  is  of  the 
same  order  of  magnitude  as  the  Da  number,  the  addition  of  diffusion  can  alter  channel 
spacing  and  lower  the  growth  rate  but  cannot  inhibit  the  instability  from  growing. 
When  diffusion  is  weak  (Pe  >  Pa),  it  can  be  neglected  altogether. 


3.6  Discussion 

The  results  of  this  work  show  the  potential  for  an  existence  of  a  channeling  instability 
in  Earth’s  upper  mantle.  The  instability  stems  from  combined  chemical  and  hydro- 
dynamical  effects.  Melt  decompresses  as  it  ascends  through  the  mantle,  increasing  its 
ability  to  dissolve  the  surrounding  matrix.  A  small  perturbation  in  porosity  allows 
a  Icirger  volume  of  unsaturated  fluid  to  flow,  thus  increasing  dissolution  and  further 
increasing  porosity  in  a  positive  feedback  mechanism. 

We  consider  a  model  system  incorporating  porous  flow,  dissolution,  and  matrix 
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compaction  eiFects  and  find  a  steady  state  formed  by  the  competition  between  disso¬ 
lution  and  compaction.  Linear  stability  analysis  is  then  performed,  predicting  that 
IoJ^Sj  113ITOW  channels  may  spontaneously  form.  The  emerging  horizontal  wavelength 
Ax  is  shown  to  depend  on  the  chemical  equilibration  length  Leq  and  the  vertical  ex¬ 
tent  of  the  system  L.  For  chemical  equilibration  lengths  that  are  smaller  than  the  size 
of  the  system  (Da  >  1),  A*  =  2.534.^ZeqL.  For  a  system  of  the  size  of  the  melting 
region,  Z-  «  75  km,  channels  will  range  from  10  cm  wide,  when  Leq  10  m,  to 
nearly  a  kilometer  wide  if  leq  ~  10  m.  Our  linear  prediction  of  channel  width  is 
in  agreement  with  (1)  postulated  vein  sizes  needed  to  successfully  segregate  the  melt 
in  order  to  see  observed  chemical  signatures  [79],  and  (2)  the  size  of  dunite  dissolu¬ 
tion  features  found  in  field  observations  [10,  19,  62,  55,  34].  This  further  strengthens 
our  belief  that  the  channeling  instability  plays  a  crucial  role  in  focusing  of  melt  and 
determining  the  geochemical  composition  of  upwelling  liquids  and  residual  mantle 
peridotites.  Channels  formed  by  the  reactive  infiltration  instability  may  provide  the 
means  for  extracting  melt  out  of  the  viscously  deformable  upper  mantle. 

If  the  chemical  equilibration  length  is  too  long,  comparable  to  the  compaction 
length,  formation  of  the  instability  is  inhibited  by  compaction.  A  coexisting  solution 
predicts  unstable  growth  of  elongated  traveling  porosity  and  concentration  waves 
that  exist  under  all  conditions.  The  exact  criterion  for  stationary  channels  to  grow 
is  Cfc^/Ireq  >  1-  We  test  whether  this  condition  is  met  in  Earth’s  mantle,  using  a 
range  of  numbers  explained  in  Appendix  B  and  tabulated  in  Table  3.1.  Equilibration 
lengths  range  from  <1  mm  to  10  m,  and  we  use  =  2  x  10"®  m"^;  c^/i/Leq  ranges 
from  0(1)  to  0(10®)  when  the  compaction  length  is  taken  as  1000  m  and  from  0(10-^) 
to  O(IO^)  when  a  compaction  length  is  taken  as  100  m.  Parameter  values  that  are 
believed  to  be  characteristic  of  Earth’s  mantle  are  thus  mostly  in  the  regime  that 
eJlows  for  rapid  growth  of  stationary  channels. 

The  effect  of  diffusion  is  also  investigated.  W^hen  diffusion  is  strong  enough,  it  will 
cause  an  increase  in  channel  spacing  and  a  decrease  in  channel  growth  rate.  However, 
this  will  not  inhibit  the  channels  from  growing. 

The  forming  channels  are  elongated  and  span  the  vertical  box  size.  Their  ampli- 
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tude  increcises  with  height,  thus  showing  increasing  chemical  disequilibration  as  liquid 
ascends  in  the  melting  column.  When  Da  is  high,  chemical  disequilibration  is  small, 
but  dissolution  of  channel  is  most  efficient.  Hence  soluble  phases  may  become  reduced 
in  volume  or  completely  exhausted.  When  channels  thus  become  stripped  of  soluble 
phases  (i.e.,  stripped  of  pyroxene  but  not  of  olivine),  melt  flowing  through  them  will 
be  effectively  isolated  from  equilibration  with  these  phases.  Since  the  channels  span 
the  system  vertically,  they  can  bring  to  the  surface  melt  that  has  not  equilibrated 
with  its  surroundings  since  it  began  its  ascent  at  the  bottom  of  the  melting  column. 
This  may  explain  why  MORE  is  out  of  chemical  equilibrium  with  upper  mautle  peri- 
dotites  [28]  and  includes  chemical  signatures  from  the  bottom  of  the  melting  column 
[37,  68]. 

Channels  formed  as  a  result  of  this  instability  should  be  identifiable  in  the  geo¬ 
logic  record.  Their  contact  relationships  should  be  indicative  of  replacement  of  host 
peridotite  as  a  result  of  selective  dissolution  of  more  soluble  phases  (i.e.,  pyroxene). 
Additionadly,  minerals  in  dissolution  chcinnels  should  be  close  to  equilibrium  with 
MORB  and  therefore  will  have  very  different  minor  and  trace  element  compositions 
from  the  same  minerals  in  surrounding  peridotite.  In  fact,  these  characteristics  are 
observed  in  dunite  bodies  within  the  residual  mantle  peridotite  section  of  the  Oman 
ophiolite  [32]. 

However,  we  emphasize  that  the  results  of  the  present  study  may  not  be  directly 
applicable  to  the  mantle.  In  particular,  the  morphology  of  dissolution  channels  arising 
from  this  instability,  developing  over  finite  length  scales  and  timescales,  cannot  be 
predicted  from  linear  stability  analysis.  It  is  possible  that  a  much  richer  behavior 
emerges  due  to  non-linear  effects.  For  example,  space  or  time  dependent  behavior,  as 
tentatively  sketched  in  Figure  3-9.  Additionally,  the  multicomponent  melt  migration 
process  in  the  earth  with  the  inclusion  of  advective  heat  transport  and  background 
melting  effects  is  more  complicated  than  the  problem  studied  here  and  is  a  topic  for 
ongoing  studies. 

In  addition  to  the  formation  of  stationary  dissolution  channels,  another  coexist¬ 
ing  linear  solution  indicates  the  existence  of  traveling  compaction-dissolution  waves. 
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whose  amplitude  increases  with  time.  These  unstable  waves,  pockets  of  undersatu¬ 
rated  melt  in  a  high  porosity  region,  will  transport  increasingly  undersaturated  melts 
to  the  surface  even  when  the  matrix  is  too  weak  to  support  stationary  chcinnels.  They 
may  also  aid  in  explaining  the  production  and  growth  of  “magmons,”  compaction 
waves  previously  thought  to  be  generated  by  an  initial  step  in  porosity. 
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Vciriable  Symbol  Vcdue _ Notes 


Solubility  gradient 

P 

2  X  10-®m-^ 

1 

Linear  dissolution  rate 

10-”  -  10-®  m  s-i 

2 

Solid  density 

P» 

3  X  10^  kg  m"^ 

Reaction  rate  constant 

R 

10“®  —  10“®  kg  m“^  s 

-1 

Grain  edge  length 

d 

lO-'*  -  10-^  m 

Total  surface  area 

10^  -  10®  mVm® 

3 

Porosity 

<f> 

10-3  -  10-2 

4 

Solid/liquid  surface  area 

S 

(1  -3)  X  10-^ 

5 

Volume  fraction  of  soluble  phase 

10-2  -  10-^ 

6 

Permeability  exponent 

n 

2-3 

7 

Melt  fraction 

F 

0.05  -  0.2 

8 

Solid  upwelling  rate 

Vo 

10-10  _  iq-9  g-l 

9 

Background  fluid  velocity 

Wo 

10-®  -  10"®  m  s-^ 

9 

Equilibration  length 

Leq 

10-^  -  10  m 

9 

Damkohler  number  (L=100m) 

Da 

10®  -  1 

10 

Diffusion  coefficient 

D 

10-12  _  10-10  m2  s-^ 

11 

Peclet  number  (L=100m) 

Pe 

10®  -  10® 

12 

compaction  length 

h 

100  —  1000  m 

13 

Table  3.1:  Characteristic  Values  Believed  to  Be  Applicable  to  Earth’s 

Mantle. 

detailed  discussion,  see  Appendix  B.  Notes:  l.[36,  35].  2.  higher  dissolution  rates  are 
from  [12]  and  [41],  and  low  dissolution  rate  is  from  [96].  3.  using  d,  cubic  grains,  and 
[90].  4.  [27,  75].  5.  in  fraction  of  total  surface  area.  6.  i.e.,  fraction  of  pyroxene  [34]. 
7.  [90].  8.  fraction  of  solid  mass  that  has  melted.  9.  see  Appendix  B.  10.  using  range 
of  ieq  given  above.  11. [25].  12.  using  wo,D  given  above.  13.  [44]. 
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Figure  3-1:  Setup  of  the  problem:  Fluid  is  driven  upward  by  a  pressure  gradient, 
entering  the  soluble  porous  media  at  z  =  0  and  leaving  at  z  =  L.  Owing  to  decom¬ 
pression,  the  fluid  has  increasing  ability  to  dissolve  the  solid  matrix.  The  matrix  is 
allowed  to  contract  by  compaction  and  thus  counteract  to  some  degree  the  effects  of 
increasing  dissolution. 
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Figure  3-2:  Nondimensional  steady  state  concentration  as  a  function  of  height,  drawn 
as  a  solid  line.  After  a  narrow  boundary  layer  («  1  reaction  length),  the  deviation  of 
the  steady  state  concentration  from  the  equilibrium  value  (represented  by  a  dashed 
line),  approaches  a  constant  value  Ac  ~  IjDa. 
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Figure  3-3:  Nondimensional  growth  rate  a  versus  nondimensional  horizontal 
wavenumber  I,  for  Da  =  10  in  a  rigid  medium  (1/a  =  0),  and  in  a  compacting 
medium  with  a  =  1;  Imeix  indicates  the  fastest  growing  wavenumber.  Comparing 
the  growth  of  unstable  channels  in  a  compacting  and  a  rigid  medium,  one  notes  that 
the  most  unstable  wavelength  is  hcirdly  ciltered  and  its  growth  rate  is  only  slightly 
lowered  due  to  the  stabilizing  effect  of  compaction.  Compaction  does,  however,  damp 
the  long-wavelength  perturbations.  This  effect  leads  to  a  critical  Da  for  existence  of 
the  instability  in  a  compacting  medium. 
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Figure  3-4:  (a)  The  fastest  growing  horizontal  wavelength,  Ax  =  27r//niax,  derived 
from  plots  similar  to  Figure  3-3  with  different  Da  numbers  in  a  rigid  medium.  For 
Da  >  1,  and  spanning  4  orders  of  magnitude,  A^  cx  If  Da.  (b)  Aspect  ratio  of 
channels,  qjl  —  Ax/A*,  cis  a  function  of  Da  in  a  rigid  medium.  Perturbations  become 
more  finger  like  as  Da  number  is  increased,  (c)  Growth  rate  of  /max  as  a  function 
of  Da  in  a  rigid  medium.  The  system  becomes  increasingly  unstable  with  increase  in 
Da  and  reaches  a  constant  limit  for  Da  — >  oo.  This  limit  is  set  by  the  derivative  of 
permeability  with  respect  to  porosity. 
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Figure  3-5:  Eigenfunctions  for  a  compacting  medium,  Z?a  =  10,  a  =  1  (a)  Pertur¬ 
bation  pressure  as  a  function  of  x  and  z.  Dark  shading  indicates  negative  valuesj 
light  shadings  indicate  positive  values.  Owing  to  the  boundary  conditions,  pressure 
perturbations  are  forced  to  zero  at  the  top  of  the  box,  thus  attaining  a  maximum 
amplitude  at  z  ~  3/4.  (b)  Porosity,  vertical  velocity,  and  concentration,  which  look 
very  similar,  are  plotted  as  a  function  of  x  and  z.  These  variables  have  increasing 
amplitude  as  a  function  of  z  and  attain  a  maximum  at  z  =  1. 


Figure  3-6:  (a)  Schematic  drawing  of  a  planar  wave,  with  wave  vector  =  P  +  w?, 
oriented  at  angle  6  to  the  vertical.  The  phase  velocity  is  the  velocity  at  which  the 
pha.ses  propagate,  (b)  (top)  Growth  rate  cr-ji  and  phase  velocity  Cp  of  compaction- 
dissolution  waves  as  a  function  of  the  planar  wave  orientation  9.  Calculations  were 
made  using  K  =  10,  Da  =  10,  and  a  =  1.  Phase  velocity  is  maximum  when  the  wave 
is  at  some  angle  to  the  vertical  and  is  zero  for  vertical  stripes.  Although  these  stripes 
do  not  propagate,  they  grow  the  most  rapidly,  as  seen  from  (t(9).  The  stationary 
vertical  stripes  are  the  collapse  of  the  plane  wave  solution  to  the  unstable  channels 
solution,  (bottom)  The  ctj  and  Cp  as  a  function  of  Da.  Calculations  were  made 
using  if  =  10,  0  =  7r/4  and  a  =  1.  The  maximum  in  phase  velocity  is  for  waves 
with  K  cos  9  =  Da.  At  this  point  a  rapid  transition  from  nearly  zero  growth  rate  to 
maximum  growth  rate  can  be  observed. 
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Figure  3-7:  Growth  rate  a  as  function  of  horizontal  wavenumber  /  for  a  rigid  matrix 
with  Da  =  10  for  various  diffusion  coefficients.  Solid  line  is  the  growth  rate  in  the 
case  of  no  diffusion,  replotted  from  Figure  3-3.  Solid  circles  are  the  linear  stability- 
analysis  results  for  a  rigid  system  with  weak  diffusion,  Pe  =  100.  Strong  diffusion 
with  Pe  =  10  is  shown  in  crosses.  The  horizontal  wavelength  (=  27r/lniax)  of  forming 
channels  becomes  larger  in  the  presence  of  diffusion  and  the  growth  rate  is  somewhat 
lowered,  but  the  nature  and  the  robustness  of  the  instability  are  not  changed.  These 
results  are  in  agreement  with  scaling  arguments. 
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Figure  3-8:  The  fastest  growing  horizontal  wavelength,  A*  =  27r//inax,  derived  from 
plots  similar  to  Figure  3-7  with  different  Pe  numbers,  for  a  rigid  medium  with  a 
constant  Da  =  10,  plotted  by  solid  circles.  In  the  case  of  strong  diffusion,  when 
Pe  <  Da,  A^  oc  1  j Pe.  On  the  same  plot  we  show  the  growth  rate  a  by  open  squares, 
to  be  unciffected  by  diffusion  when  Pe  is  large  but  to  be  reduced  when  diffusion 
becomes  important.  Channels  will  grow  {a  >  0)  for  all  diffusion  rates. 
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Figure  3-9:  When  high-porosity  channels  form,  the  solid-liquid  surface  area  per  unit 
volume  is  reduced  and  the  characteristic  equilibration  length  of  the  system  will  in¬ 
crease.  This  means  that  the  wavelength  of  the  channels  will  increase  as  well  (from 
equation  (55)).  If  we  also  consider  the  increase  in  melt  velocity  with  height,  we  get 
an  additional  increase  in  Teq.  Thus  we  can  tentatively  propose  the  following  picture: 
narrow  channels  that  form  deep  in  the  melting  column  will  act  as  the  background 
porous  structure  for  a  higher  level  and,  due  to  the  increase  in  horizontal  wavelength 
with  equilibration  length,  will  coalesce  to  form  an  interconnected  network  of  high 
porosity  tubes  similar  to  the  upside  down  “fractal  tree”  proposed  by  [22]. 
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Chapter  4 


Simulations  of  flow  and  reaction 
in  porous  media 


Abstract 

I  present  a  3D  computer  model  that  simulates  flow  and  reaction  through  a  porous 
medium,  by  solving  the  relevant  mass  conservation  partial  differential  equations  and 
Darcy’s  law.  Reaction  generally  results  in  re-organization  of  the  porous  media.  Dis¬ 
solution  causes  focused  and  highly  correlated  flow  patterns,  while  deposition  reduces 
correlations  in  the  spatial  cirrangement  of  porosity,  so  as  to  cause  dispersion  of  flow. 
Three  specific  experiments  relevant  to  melt  migration  in  the  mantle  are  also  per¬ 
formed.  The  first  simulates  conditions  in  upwelling  melt  beneath  mid-ocean  ridges, 
where  I  find  formation  of  high  permeability  channels  cis  predicted  in  Chapter  III.  The 
second  experiment  simulates  conditions  for  upwelling  melt  beneath  intra-plate  vol¬ 
canos  or  magmatic  arcs,  where  melt  crystalization  is  expected.  Simulations  of  such 
conditions  indicate  increasingly  diffuse  and  homogeneous  porous  flow,  even  when  flow 
was  initially  channelized.  The  third  experiment  is  of  melt  changing  from  corroding 
the  matrix  to  depositing  in  the  matrix.  In  the  transition  zone  a  low  porosity  cap 
forms  over  a  high  porosity  region,  with  rapidly  increasing  overpressurization.  The  re¬ 
sults  have  implications  for  degree  of  chemical  equilibration  of  melt  in  different  regions 
as  well  as  for  domincint  modes  for  melt  migration  (i.e.,  whether  porous  channels  or 
periodically  forming  cracks  will  appear). 


4.1  Introduction 

In  this  paper  I  present  a  newly  developed  3D  numerical  model  of  flow  and  reaction 
through  a  porous  media,  and  initial  results  from  it.  There  are  not  many  studies  about 
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coupled  flow  and  reaction,  even  less  about  deposition,  and  none  of  the  existing  ones 
handles  both  deposition  and  dissolution.  This  work  aims,  first  of  all,  to  study  the 
evolving  nature  and  statistics  of  reacting  porous  media  and  to  establish  the  differ¬ 
ences  between  dissolution  and  deposition  as  opposite  processes.  The  secondary  aim 
of  this  work  is  to  perform  some  initicii  numerical  experiments  mimicking  different 
environments  for  upwelHng  melt  in  the  mantle,  and  compare  the  results  to  existing 
observations  and  predictions. 

The  method  of  study  involves  solving  a  set  of  partial  differential  equations,  which 
describe  mass  conservation  of  fluid  and  rock,  and  Darcy’s  law.  This  mathematical 
description  is  applicable  to  scales  above  the  pore  size,  where  one  can  define  an  average 
porosity  and  permeability.  The  resulting  set  of  coupled  equations  is  solved  on  a 
computer  (a  SUN  workstation)  using  a  Multi-grid  method  [61]  to  obtain  the  pressure 
field.  This  efficient  method  enables  inclusion  of  up  to  64^  nodes  on  a  workstation, 
where  memory  availabihty  limits  simulations  of  larger  boxes. 

Before  deciding  upon  a  macroscopic  partial-differential-equation  solver  I  have  at¬ 
tempted  two  additional  numerical  methods,  which  I  will  spend  a  brief  paragraph  to 
describe,  in  the  dim  chance  that  it  will  save  somebody  some  work  in  the  distant  fu¬ 
ture.  First  I  constructed  a  version  of  the  lattice-gas  method  [66]  in  which  particles 
flow  and  react  on  a  lattice.  This  method  proved  to  be  well  suited  for  microscopic  sim¬ 
ulations.  The  initial  results  were  easy  to  obtain  and  provided  qualitative  indicators, 
but  I  was  unable  to  continue  into  quantitative,  large  scale  studies,  due  to  size  limi¬ 
tations.  The  description  of  the  model  and  the  initial  results  are  presented  in  [36].  I 
then  proceeded  to  implement  a  network  model  in  which  a  reactive  fluid  flows  through 
pipes  of  various  sizes  and  may  change  the  radius  of  the  pipes  by  chemical  reaction. 
At  first  gljmce  this  method  seemed  most  promising,  but  it  is  limited  by  the  need  to 
solve  for  transport  of  solute.  This  is  done  either  by  particle  tracking  methods,  which 
limit  the  size  of  the  system  to  be  investigated  similar  to  the  Lattice- Gas  and  other 
microscopic  methods,  or  by  differential  equation  solutions,  which  have  the  limitations 
of  macroscopic  models.  Thus  I  concluded  that  for  quantitative  study  of  large  scale 
organization  a  macroscopic  model  that  solves  a  set  of  partial  differential  equations 
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would  provide  the  most  useful  tool. 

In  the  first  section  of  this  chapter  I  mathematically  describe  reactive  fiow  through 
a  porous  medium,  similar  to  the  macroscopic  description  in  Chapter  III  [3],  but  here  in 
the  rigid  medium  limit.  In  the  second  section  the  computer  method  is  described  and 
validated.  In  the  final  section  I  present  two  different  sets  of  experiments  of  dissolution 
and  deposition  in  a  porous  media. 

The  first  set  of  experiments  constitutes  a  basic  study  of  3D  flow  and  reaction 
in  a  porous  media.  The  statistical  and  general  properties  of  the  porous  media  are 
investigated,  after  a  dissolving  or  depositing  fluid  have  flown  through  it.  Results  from 
dissolution  and  deposition  axe  compcired  to  one  another  and  to  the  initially  randomly 
constructed  media. 

The  second  set  of  experiments  is  aimed  at  mimicking  (very  simplisticly)  conditions 
of  melt  ascending  beneath  (1)  mid-ocean  ridges,  similar  to  the  situation  described 
in  Chapter  III,  and  (2)  subduction-related  magmatic  arcs  or  intraplate  hot  spots. 
Geochemical  observations  on  lavas  and  mantle  samples  indicate  that  melt  ascending 
beneath  mid-ocean  ridges  is  out  of  chemical  equilibrium  with  surrounding  mantle  [28, 
27],  while  observations  of  lavas  and  mantle  samples  from  hot  spots  and  arcs  genereilly 
indicate  extensive  reactions  (for  review  see  [36]).  Such  geochemical  observations  are 
accompcinied  by  geologiccil  observations  of  different  structures  for  melt  transport  in 
the  two  cases:  replacive  dunites  in  the  case  of  mid-ocean  ridges  and  fractures  for  intra¬ 
plate  basalt.  This  set  of  experiments  is  a  continuation  of  work  initiated  in  [36]  and 
Chapter  III,  in  which  the  chemical  and  geological  observations  are  explained  by  the 
different  chemical  interactions  that  occur  between  melt  and  the  surrounding  matrix: 
In  midocean  ridges,  dissolution  channels  form  when  melt  upwells  in  close  to  adiabatic 
conditions  [36,  3,  32].  In  hot  spots  and  arcs,  melt  must  pass  through  a  conductively 
cooled  tectosphere  and  so  crystallize,  possibly  resulting  in  an  overpressurized  and 
finally  cracked  media  [36].  In  addition  to  checking  the  flow  patterns  in  these  two 
different  Ccises,  I  also  simulate  the  transition  between  dissolution  and  deposition, 
(a  trcinsition  likely  to  happen  as  melt  passes  from  a  an  adiabatic  to  a  conductive 
geotherm)  and  discuss  the  likelihood  and  position  of  fracture  initiation. 
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All  simulations  were  performed  with  relatively  high  reaction  rates  and  low  diffu- 
sivities  so  as  to  maximize  observed  trends  and  minimize  simulation  times.  The  sets 
of  experiments  are  by  no  means  exhaustive.  For  discussion  of  the  effect  of  varying 
the  different  control  parameters  see  Chapter  III,  [80,  67]  and  [9]. 


4.2  Formulation  of  the  problem 

The  description  derived  below  is  close  to  previously  used  forms  [14,  57,  3],  but  here 
an  analytical  description  of  deposition  is  added.  The  mathematical  formalism  uses  a 
conventional  continuum  approach  suited  for  description  of  phenomena  that  occur  at 
and  above  the  Darcy  scale. 

4.2.1  General  Equations 

Mass  conservation  equations  are: 


dp.il  -  <f>) 
dt 

=  Er. 

i 

(4.1) 

+  V  •  (P/v^) 

-  -Er.. 

(4.2) 

dp,il  -  4>)c\ 
dt 

=  v-[D,v,(i-,^)Vc']  +  ri, 

(4.3) 

+  V  •  (p/^VC.) 

=  V  •  [DiPf(f)Vci]  -  r,-. 

(4.4) 

where  (f>  is  the  porosity,  p,  and  pf  are  the  solid  and  fluid  density  respectively  in  kg 
m~^,  and  T,-  is  the  mass  transfer  rate  of  mineral  i  from  fluid  to  solid  in  kg  m“^  s“  1.  D,- 
and  are  the  diflusion  coefflcients  of  component  i  in  the  fluid  and  solid  respectively 
snd  CijC^  are  the  mass  fractions  of  component  t  in  the  fluid  and  solid  respectively, 
with  'Zi  c.-  =  1  and  'Zi  c-  =  1. 

Equations  4.1  and  4.2  describe  the  conservation  of  total  solid  and  fluid  mass  re¬ 
spectively.  Equations  4.3  and  4.4  describe  the  conservation  of  each  mineral  component 
in  the  solid  and  fluid  phase  respectively. 

Assuming  first-order  chemical  reaction  kinetics,  the  rate  of  mass  transfer  of  com- 
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ponent  i  between  fluid  and  solid  is 

Fi  =  RiAi{x, t){ci  -  c^q.),  (4.5) 

where  Ri  is  the  reaction  rate  constant  of  component  i  in  kg  s“^  and  Ai(x,  t)  is 
the  effective  speciflc  surface  area  (m’/m^)  available  for  reaction  of  component  i- 
is  the  equilibrium  concentration  of  minertil  i  in  the  fluid,  given  in  mass  fraction.  Cg,^ 
may  be  a  function  of  position  and  time,  as  in  the  case  for  varying  temperature  or  pH 
conditions.  The  specific  surface  area,  Ai(x,  t)  is  also,  in  general,  a  complex  function 
which  has  different  values  depending  on  whether  a  mineral  is  dissolving  or  depositing, 
and  may  depend  on  the  presence  of  other  minerals.  Here  I  shall  approximate  Ai{x,  t) 
in  a  simple  form.  The  effective  specific  surface  area  for  dissolution,  termed  Ai^“*, 
is  a  function  of  the  amount  of  mineral  i  available  for  dissolution  at  the  pore-grain 
interface.  Since  M*i  =  -  (f))  is  the  mass  fraction  of  solid  component  i, 

Af«(c',  4>)  oc  oc  (c.i(l  -  <f>)Y^\  (4.6) 

Equation  4.6  arguably  applies  to  dissolution  and  not  to  deposition,  since  dissolution  is 
limited  by  the  mass  of  soluble  material  in  the  solid  phase,  while  deposition  is  generally 
independent  of  the  amount  of  mineral  already  existing  in  grains,  but  depends  only 
on  the  total  surface  area  exposed  in  the  pores.  Since  the  volume  fraction  occupied 
by  pores  is  <f>,  the  specific  surface  area  associated  with  this  volume  fraction  will  be 
approximately 

oc  (4.7) 

The  closing  equation  is  Darcy’s  law,  which  relates  the  pressure  p  to  the  flux  of  the 
fluid: 

(pv  =  — — Vp,  (4-8) 

where  p  is  the  viscosity  of  the  fluid,  p  is  the  pressure  in  excess  of  hydrodynamic 
pressure,  and  k  is  the  permeability. 

The  relation  between  the  permeability  and  the  geometry  of  the  matrix  is  important 
for  correct  modeling  of  the  non-linear  process  involved  in  this  study,  but  a  clear 
analytical  description  is  not  obvious.  For  example,  it  is  possible  that  by  selective 
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deposition  (or  dissolution)  in  strategically  positioned  necks  the  porosity  will  be  hardly 
changed,  while  the  permeability  will  be  highly  altered.  Never  the  less,  in  the  literature 
(e.g.  [8,  88])  permeability  is  usually  tciken  to  be  a  power  law  function  of  the  porosity 

A:(^)  =  (4.9) 

with  d  a  typical  grain  size,  n  usually  between  2  and  3,  and  h  a  constant.  This  relation 
is  not  very  good  for  low  porosities,  but  is  used  in  this  study  for  lack  of  a  better 
description,  and  as  a  first,  even  though  not  perfect,  step.  This  approximation  can  be 
peirtiaJly  justified  when  one  realizes  that  the  discretized  version  of  the  set  of  equations 
presented  above  (4. 1-4.9)  is  identical,  when  n  =  2,  to  the  set  needed  for  solving  for  flow 
and  reaction  through  a  network  of  tubes.  Thus,  similarly  to  a  network  model,  each 
node  (in  the  network  model  representing  a  tube)  by  itself  obeys  4.9,  but  the  global 
arrangement  in  space  may  result  in  macroscopic  porosity-permeability  relations  that 
are  different  than  4.9. 

Boundary  conditions.  In  general,  equations  4. 1-4.9  will  require  either  pressure  or 
flux  boundary  conditions  on  both  the  total  fluid  mass  and  on  each  individual  fluid 
component.  Initial  conditions  on  <}>  and  Ci  are  also  needed  in  order  to  determine  the 
spatial  and  temporal  evolution  of  the  six  (in  the  case  of  a  single  soluble  component 
in  ZD  flow)  unknowns  y,  p,  (l>  and  c. 

4.2.2  Simplifications 

For  simplicity,  I  eissume  that  the  solid  phase  is  composed  of  two  components:  a  soluble 
component  which  can  chemically  react  with  the  fluid  by  dissolution  or  precipitation, 
and  a  insoluble  component.  The  soluble  material  has  mass  fraction  c,  in  solid  and 
reaction  rate  R  while  the  non-soluble  component  has  mass  fraction  1  —  c,  in  solid  and 
its  reaction  rate  is  0.  Diffusion  within  the  solid  will  be  considered  negligible.  The 
fluid  phase  is  composed  of  a  Ccirrier  fluid  with  mass  fraction  1  -  c  and  a  dissolved 
component  with  mass  fraction  c.  The  carrier  fluid  does  not  enter  the  solid  phase. 
The  density  of  the  fluid  and  solid  phases  are  presumed  constant  as  the  composition 


102 


of  the  melt  changes  due  to  chemical  reaction.  The  subscript  i  will  be  dropped  from 
this  point  on. 

In  this  simplified  case  we  can  use  4.3  to  write  the  mass  conservation  equation  of 
the  non- soluble  phase, 

(1  —  C4)(l  —  <f>)  =  Vn,  =  const.  (4-10) 


From  4.10,  4.6  and  4.7  we  write  the  effective  surface  area  as  a  function  of  the  porosity 
only: 


A{4>)  OC  < 

1  (^/  -  if  C  <  Cg, 

[  if  c>  Cg, 

(4.11) 

where  (f>f  =  1  —  Vnt  would  be  the  final  porosity  if  all  soluble  material  is 
Equations  4. 1-4.4  can  now  be  written  as 

dissolved. 

d<f> 

dt 

(4.12) 

|  +  V.(v« 

=  -r//>/. 

(4.13) 

,dc 

4>^  +  •  Vc 

=  DV-{<f>Vc)-{l-c)Vlpf, 

(4.14) 

r 

=  RA{(l>){c- Ceq), 

(4.15) 

where  4.14  is  a  result  of  subtracting  4.2  from  4.4  and  equation  4.3  has  been  used  to  de¬ 
termine  the  reactive  surface  area  for  dissolution,  and  will  not  be  followed  individually 
anymore. 

4.2.3  Nondimensionalization. 

In  nondimensionalizing  equations  4.11-4.15  and  4.8,  I  define  characteristic  effective 
permeability  as  A:o  =  characteristic  flux  as  <f>oWQ  =  kopof  fiL,  where  L 

is  a  length  scale  of  interest,  and  po  is  the  characteristic  pressure  difference  across  L. 
A  chciracteristic  concentration  is  defined  as  cq  =  Cg,,  which  scales  the  concentration 
c  to  the  spatial  and  temporal  average  of  the  equilibrium  concentration  Cg,.  This 
definition  reduces  to  Co  =  Cg,  when  the  equilibrium  concentration  is  constant.  Finally 
a  chciracteristic  surface  area  for  reaction  will  be  defined  as  Aq  =  . 
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I  next  define  the  nondimensional  variables,  denoted  by  primed  letters: 


X  =  Lx' 

<t>  = 

P  =  Pop' 

V  =  Wqv' 


Wq 

C  =  Cqc' 

Ce?  =  Co/(x,i) 

<f>f  =  4>0<f>'f 


(4.16) 

(4.17) 


Equation  4.16  defines  the  Damkholer  number,  which  measures  the  advection  time 
scale  versus  the  reaction  time  scale  in  the  system.  Equation  4.17  defines  the  Peclet 
number,  which  mecisures  diffusion  time  scales  versus  the  advection  time  scales  in  the 
problem. 

Using  the  non- dimensioned  variables  and  parameters,  4.8,  4.9  and  4.11-4.15  can 
be  rewritten,  eifter  some  algebra,  as: 

d<f>'  Pf 

SF  = 

V.(rvp')  =  c„^^.^Dar'  (4.19) 

dd  1 

=  ;^V.(^'Vc')-(l-coc')Dar'  (4.20) 

r  =  A{<i>'){c' -  f{-x.,t))  (4.21) 


K4>)  Oil  \  '  4.22) 

\cl>'y^  \id>f{x,t) 

This  constitutes  the  basic  set  of  equations  for  this  study.  From  this  point  on  primes 
will  dropped  with  reference  to  the  non-dimensional  variables.  Equation  4.18  describes 
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the  temporal  evolution  of  porosity  due  to  reaction.  Equation  4.19  describes  total  fluid 
conservation  which  is  simply  the  Laplace  equation  for  the  pressure.  Equation  4.20 
states  that  mineral  concentration  changes  due  to  reaction,  diffusion  and  advection. 
Finally  Equation  4.21  and  4.22  are  used  to  evaluate  the  mass  transfer  rate. 


4.3  Description  of  numerical  model 

The  numerical  method  for  solution  of  4.18-4.22  uses  an  explicit  partial  differential 
equation  solver  which  is  implemented  on  an  evenly  spaced  grid  with  grid-spacing  A 
in  each  of  the  respective  dimensions,  and  time  steps  of  size  At.  The  model  uses  wrap¬ 
around  boundary  conditions  in  the  x  and  y  directions,  and  allows  for  flux  or  pressure 
specifications  in  the  z  =  0,L  boundaries.  The  model  is  implemented  as  follows; 

1.  Given  ZD  concentration  and  porosity  fields,  find  the  mass  transfer  rate  F  ev¬ 
erywhere  using  4.21  and  4.22. 

2.  Given  pressure  or  fliix  boundary  conditions  and  F,  find  the  ZD  pressure  field 
from  equation  4.19.  The  solution  of  such  Laplace  equations  are  most  time  and 
memory  consuming,  and  constitute  the  step  that  limits  the  size  of  the  problem 
to  be  solved.  The  model  uses  the  recently  developed  Multi-Grid  Method  [61], 
which  requires  a  number  of  operations  of  order  n,  where  n  is  the  number  of 
grid  points.  This  method  is  able  to  solve  larger  systems  than  direct  solvers  that 
need  a  number  of  operations  of  order  n^,  but  unlike  the  direct  solvers  it  cannot 
resolve  high  frequency  features. 

3.  Operator  splitting  [61]  is  then  used  to  solve  for  the  new  solute  concentration, 
from  equation  4.20.  The  advective  contribution  is  calculated  using  the  “cor¬ 
rected  up-wind  scheme”  [74],  a  stable  advection  scheme  with  courant  condition: 
AfjQl/A  <  l/\/3  where  Q  =  <^"Vp.  The  corrected  upwind  scheme  has  small 
implicit  diffusion  ptirallel  to  the  flow  direction.  The  source  (or  sink)  term  due 
to  chemical  reaction  is  given  by  F,  calculated  in  step  1.  The  diffusion  step  is 
implemented  using  a  first-order  explicit  space-centered  scheme  [61]  (stable  for 
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6Ai/ A  <  Pe).  The  stability  requirement  is  easily  met  for  the  small  diffusivities 
(large  Pe)  which  I  use. 

4.  Finally  the  new  porosity  field  is  obtained  from  4.18,  and  the  loop  returns  to 
step  1. 

In  the  validation  tests  and  the  experiments  described  below,  fluid  is  driven  by  a 
pressure  gradient  from  z  =  0  to  z  =  L,  with  constant  pressure  boundary  conditions. 


4.4  Testing  of  numerical  model 

In  order  to  test  the  model  I  had  to  concentrate  on  simplified  cases  with  many  symme¬ 
tries,  due  to  the  lack  of  theoretical  work  in  more  complex  cases.  First,  each  feature 
of  the  model  was  checked  by  itself.  Advection  and  diffusion  of  an  initial  Gaussian 
were  tested  to  give  the  analytically  predictable  profile.  This  test  did  not  incorporate 
reaction.  Reaction  in  a  uniform  porous  media  was  tested  against  the  analytically  pre¬ 
dicted  decay  to  equilibrium  of  the  concentration  field.  This  test  did  not  incorporate 
advection. 

For  a  dissolving  porous  media  there  are  a  few  cases  which  incorporate  both  reaction 
and  transport,  cind  for  which  there  are  anadytical  predictions.  I  have  tested  my  model 
against  a  few  of  these  cases.  The  test  I  present  is  a  comparison  of  the  numerical 
and  predicted  dispersion  curves  for  the  reactive  infiltration  instability.  The  analytical 
results  are  derived  for  reaction  and  flow  in  a  two  dimensional  porous  media,  and 
therefore  this  test  constitutes  the  most  demanding  of  all  available  analytical  ones. 

For  deposition  in  a  porous  media  there  are  no  available  analytical  calculations 
which  incorporate  both  reaction  and  advection,  and  so  I  will  only  present  a  test 
of  reaction  in  a  uniform  porous  media,  and  rely  on  the  fact  that  deposition  and 
dissolution  are  nearly  symmetrical  process,  and  that  the  combined  non-linear  effects 
of  the  reaction  and  transport  were  shown  to  be  consistent  with  theory  in  the  previous 
test. 
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4.4.1  Test  1:  The  reactive  infiltration  instability  in  2D 

When  an  undersaturated  fluid  is  forced  to  flow  through  a  partially  soluble  porous 
media,  a  propagating  reaction  front  will  form  [14,  57].  Upstream  of  the  reaction  front, 
where  all  the  soluble  material  has  already  been  dissolved,  the  concentration  will  be 
that  of  the  incoming  undersaturated  fluid.  Downstream  of  the  reaction  front,  the 
fluid  will  be  in  chemical  equilibrium  with  the  solid  (i.e.,  saturated)  and  no  dissolution 
will  occur,  provided  that  Da  >  1.  The  position  of  the  front  is  unstable  to  small 
perturbations:  if  the  initial  permeability  is  slightly  higher  in  a  certain  region,  this 
region  will  entrain  more  unsaturated  fluid  and  dissolve  more  rapidly,  thus  increasing 
its  porosity  even  further  in  a  positive  feedback  mechanism.  [57]  present  a  detailed 
analysis  of  this  instability  in  2D  ^  They  predict  that  for  small  perturbations,  the 
front  will  develop  scalloped  fingers  that  grow  in  an  unstable  manner  with  time,  such 
that  any  initial  sinusoidal  perturbation  in  position  will  grow  as 

C(x,  2,  t)  =  z{t)  4-  XI  cos(mi)e‘*'("*^*,  (4.23) 

m 

where  (j{m)  is  the  growth  rate  of  the  wavelength  27r/Tn  and  ^{x,z,t)  is  the  position 
of  the  front,  with  average  height  z{t).  The  resulting  dispersion  relation  is  a  function 
of  the  diffusion,  advection,  and  reaction  rates.  [57]  assume  that  Cgq  cq  constj 
cq! Pf  <<  1  cind  pf  =  Pf  Non-dimensionalization  is  performed  with  t  =  L{woCo)~^'t'- 
Their  assumptions  are  analogous  to  having  very  fast  reactions  (i.e..  Da  -*  oo),  or 
an  infinitely  thin  reaction  front.  The  analytical  prediction  for  the  linear  growth  rate 
(jj{rn)  in  this  case  is: 

“("•)  =  h + (1  - 

(^/-w(l+7)^ 

where  Vf  is  the  final  fluid  velocity  after  the  front  has  passed,  (f>b  is  the  porosity 
before  the  passage  of  the  front,  a  =  Pe*Vf  and  7  =  (f>bhl^fkf  where  h,  kj  are  the 
permeabilities  before  and  after  the  passage  of  the  front. 

To  mecisure  a;(Tn),  I  performed  a  set  of  simulations,  where  flow  enters  a  box  (33 
nodes  per  side)  at  2  =  0  with  initial  concentration  c  =  0.  An  initial  front  between  4>b 

^There  is  a  slight  mistake  in  their  paper  in  the  basic  fluid  conservation  equation,  but  since  this 
mistake  is  of  order  e  it  does  not  effect  the  final  result. 
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and  was  placed  at  ^(i,  z,0)  =  z  +  A  sin(m*s).  The  fluid  was  allowed  to  react  with 
the  porous  media  and  the  power  spectrum  of  the  porosity,  iS'(m,  t),  was  obtained.  The 
amplitude  of  the  3D  power  spectra  in  the  seeded  wavelength,  S{mx,t),  was  checked 
as  a  function  of  time.  To  derive  uj(mx)  from  such  measurements,  I  use  4.23  in  the 
following  formula:  u>{mx)  =  ln(y^(5(mx,t)/5(m*,0)))/(tco),  where  cq  is  a  parameter 
needed  to  convert  between  the  two  different  time  scales  used  in  this  paper  and  in  [57]. 
The  results  are  plotted  in  Figure  4-1. 

4.4.2  Test  2:  Deposition  in  uniform  porous  media 

Since  I  do  not  know  of  any  existing  analytical  prediction  for  combined  flow  and 
deposition  in  a  heterogeneous  porous  media,  I  performed  a  primitive  test  for  flow  and 
deposition  in  a  uniform  porous  media.  Supersaturated  fluid  fills  a  uniform  porosity 
box,  where  I  allow  for  wrap-around  boundary  conditions  in  z  as  well  as  the  usual 
wrap  around  conditions  in  x  and  y.  In  this  case  one  expects  no  spatial  concentration 
gradients.  For  coP/  <  p,,  the  change  in  porosity  at  short  time  scales  is  negligible 
(from  equation  4.18).  In  this  case  equation  4.4  has  a  solution  of 

c{t)  =  I  -  exp  ( - - - ).  (4.25) 

Figure  4-2  shows  that  simulation  results  follow  theoretical  predictions  closely. 


4.5  Simulations  and  results 

In  this  section  I  present  various  types  of  computer  experiments  and  their  results.  The 
experiments  axe  exploratory,  since  there  is  very  little  previous  quantitative  knowledge 
about  combined  reaction  and  flow  through  porous  media,  and  what  exists  is  mainly 
aimed  at  dissolution.  Many  of  the  results  will  be  quantified  by  measuring  the  correla¬ 
tion  function,  c(r),  of  different  variables  in  the  system.  The  correlation  function  in  all 
experiments  is  obtained  using  the  Weiner-Khintchine  theorem  (e.g.,  [ll]).  Practically, 
this  is  done  by  taking  the  inverse  fourier  transform  of  the  3D  power-spectra  of  the  field 
of  interest  (e.g.  porosity  or  flux).  c(r)  is  then  normalized  to  I,  and  ranges  between  —I 
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and  1.  The  correlation  function  gives  a  quantitative  indication  of  the  probability  of 
finding  identical  values  of  the  investigated  scalar  field  at  2  different  points  separated 
by  a  distcince  vector  r,  cis  explained  in  more  detail  in  the  introduction  of  Chapter  II. 

4.5.1  General  dissolution  and  deposition 

This  section  presents  a  set  of  simple  experiments.  Reactive  fluid  flows  through  a  3D 
porous  media,  which  is  wrapped  around  in  the  z  as  well  as  the  usual  x,y  directions. 
Porosity  in  the  initial  configuration  is  white  Gaussian  noise  around  a  mean.  The 
fluid  in  these  experiments  is  kept  eirtificially  from  reaching  chemical  equilibrium  in 
order  to  observe  the  effect  of  long  term  reaction.  The  fluid  is  kept  undersaturated 
(in  the  dissolution  experiments)  or  supersaturated  (in  the  deposition  experiments), 
by  readjusting  the  mean  concentration  at  each  time  step  to  a  given  constant  value. 

Figure  4-3  shows  average  porosity  within  the  box  as  a  function  of  time  for  (al) 
dissolving  porous  media  and  (bl)  clogging  porous  media.  Porosity  asymptotically 
approaches  <;6/  in  the  dissolving  case  and  0  in  the  depositing  case,  following  equation 
4.22.  Figure  4-3  also  shows  final  2D  sections  of  the  flow  field,  parallel  to  the  flow 
direction,  for  a  porous  media  subjected  to  (a2)  flow  and  dissolution  and  (b2)  flow  and 
deposition. 

Dissolution.  When  an  undersaturated  reactive  fluid  flows  through  a  rock,  one  ex¬ 
pects  focusing  of  flow  due  to  processes  that  preferentially  dissolve  regions  of  high 
initial  permeability,  similar  to  the  reactive  infiltration  instability. 

Figure  4-4  show  the  evolution  of  flux  (al)  and  porosity  (a2)  histograms  and  the 
porosity  correlation  function  (a3)  due  to  dissolution.  The  porosity  histogram  (Figure 
4-4(a2))  after  700  time  steps  of  dissolution  shows  only  a  slight  increase  in  variability 
from  its  initial  state.  On  the  other  hcind,  the  flux  histogram  (Figure  4-4(al)),  shows 
a  great  amount  of  variability  compared  to  its  initial  configuration.  This  can  be 
explciined  by  the  rearrangement  of  the  pore  distribution  in  a  correlated  way,  so  as  to 
form  elongated  high  permeability  pipes  with  high  flux,  alternating  with  regions  of  low 
permeability  and  low  flux.  The  contrasting  fast  and  slow  regions  of  flow  thus  reflect 
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the  variability  in  permeability  due  to  spatial  organization  of  porosity,  rather  than  due 
to  the  actual  change  in  porosity  itself.  This  can  be  clearly  observed  in  Figure  4-4{a3) 
as  an  increase  in  the  correlation  function  of  the  porosity  in  the  direction  parallel  to 
the  flow  (lag  in  z),  which  meeins  that  high  porosity  is  likely  to  be  found  below  or  above 
other  regions  of  high  porosity  while  low  porosity  regions  will  be  found  adjacent  to 
other  low  porosity  regions  (i.e.,  elongated  pipes).  This  increased  correlation  function 
is  most  prominent  along  the  flow  direction,  indicating  that  the  medium  has  become 
anisotropic. 


Deposition.  An  interesting  and  not  well  understood  problem  is  how  flow  reorga¬ 
nizes  as  deposition  clogs  the  porous  media.  Qualitatively,  one  expects  that  since  fluid 
flows  faster  in  regions  of  high  permeability,  it  will  tend  to  remain  further  from  equilib¬ 
rium  in  those  regions.  In  the  case  of  deposition  this  means  that  fluid  retains  a  higher 
supersaturation  in  the  high  permeability  regions,  and  thus  will  be  able  to  deposit 
more  in  these  regions.  This  process  is  expected  to  cause  a  negative  feed-back:  any 
highly  permeable  regions  are  choked  by  increased  deposition,  and  the  flow  becomes 
more  diffuse  with  time. 

This  phenomena  is  observed  in  our  simulations.  Most  of  the  deposition  occurs  in 
the  high  permeability  regions,  where  the  fluid  tends  to  clog  the  smaller  necks,  thus  ef¬ 
fectively  blocking  initially  preferred  paths  for  flow,  and  causing  flow  to  become  diffuse 
and  uniform,  as  observed  in  the  spiky  flux  histogram  in  Figure  4-4(bl).  On  the  other 
hand,  the  porosity  distribution  is  not  much  narrower  than  in  its  initial  state  (Figure 
4-4(b2)).  It  is  thus  probable  that  reduction  of  correlations  in  the  permeability  path 
are  responsible  for  flow  dispersivity  and  uniformity,  as  opposed  to  homogenization  of 
porosity.  The  correlation  function  of  the  porosity,  shown  in  Figure  4-4(b3)  demon¬ 
strates  that  correlations  along  the  direction  parallel  to  flow  are  somewhat  reduced  at 
small  z  lag. 

The  reduction  of  correlations,  and  the  possible  formation  of  negative  correlations, 
can  be  explained  in  the  following  way:  since  deposition  is  a  function  of  the  flux  of 
supersaturated  material,  it  will  occur  most  rapidly  in  the  main  paths  for  flow.  These 
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paths  are  regions  of  connected,  relatively  high  permeability,  present  in  the  initial 
random  cirrangement.  Deposition  in  these  paths  will  reduce  the  upstream  permeabil¬ 
ity  eind  preferentially  cause  regions  of  clogged  low  permeability  to  be  connected  to 
downstream  regions  of  high  permeability,  which  may  result  in  negative  correlations. 


4.5.2  Dissolution  and  deposition  in  mantle  flow 

In  this  section  I  present  results  from  simulations  directly  aimed  at  better  under¬ 
standing  flow  and  chemical  patterns  of  melt  upwelling  in  the  mcintle,  under  various 
conditions. 

The  first  experiment  is  for  fluid  increasingly  able  to  dissolve  the  porous  matrix 
through  which  it  is  flowing.  This  problem  is  set  to  mimic  aspects  of  melt  flow  in 
upwelling  mantle  beneath  mid-ocean  ridges,  where  decompression  of  ascending  mantle 
causes  an  increase  in  solubility  of  solid  phases  with  height  [72,  36].  A  2D  linear 
compacting  system  was  studied  in  Chapter  III  [3]  and  showed  the  formation  of  fast 
flowing  melt  channels.  The  3D  non-linear  problem  (in  the  rigid  limit)  is  investigated 
here. 

The  second  experiment  is  for  melt  ascending  while  crystallizing.  [36]  asserted  that 
when  melt  ascends  beneath  magmatic  arcs  and  hot  spots  volcanos  it  will  first  dissolve 
the  surrounding  mantle,  similar  to  the  situation  in  midocean  ridges.  But,  in  contrast 
to  midocean  ridges,  it  will  encounter  a  wide  region  of  conductively  cooled,  static 
mcintle  below  the  crust,  cind  so  the  melt  will  cool  to  below  the  pyroxene  saturation 
point  and  begin  to  crystallize.  Numerical  simulations  are  set  to  investigate  flow 
patterns  in  such  crystallizing  regions. 

The  third  experiment  investigates  the  transition  zone  between  dissolution  and 
crystallization.  Zones  of  transition  between  high  and  low  permeability  have  been 
proposed  for  “magma  pooling”  and  overpressurization,  and  cis  probable  places  for 
fracture  initiation  [86,  56,  36].  Kelemen  et  al.  [36]  emphasized  that  these  are  par- 
ticulcirly  likely  where  melt-rock  reaction  undergoes  a  transition  form  net  dissolution, 
under  adiabatic  conditions,  to  net  precipitation,  under  conductively  cooled  conditions. 
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Melt  fraction  increasing  upstream  (adiabatic  conditions).  The  experimental 
conditions  cire  such  that  fluid  enters  at  z  =  0  with  concentration  c  =  0.  The  enter¬ 
ing  fluid  is  in  equilibrium  with  the  surrounding  matrix.  The  equilibrium  saturation 
increases  linearly  upstream,  so  that  Co(z)  =  z/i,  where  z  is  the  height  and  L  is  the 
length  of  the  box,  (in  A).  The  loccil  concentration  is  thus  a  balcince  between  reaction 
and  advection  of  undersaturated  matericd  form  downstream,  as  explained  in  Chapter 
III.  For  Da  >•  1  the  fluid  maintains  a  close  to  equilibrium  concentration,  but  dissolu¬ 
tion  occurs  rapidly  leading  to  channel  formation  as  predicted  in  Chapter  III  for  2D. 
Figure  4-5a  shows  a  vertical  slice  (parallel  to  the  flow  direction)  of  the  ZD  flux  field  af¬ 
ter  channels  have  been  established  by  the  dissolution  process.  Figure  4-6a  and  b  show 
the  correlation  function  in  the  direction  parcillel  and  perpendicular  to  the  pressure 
gradient,  respectively.  Flow  is  focussed  and  thus  exhibits  strong  correlations  in  the 
direction  parallel  to  flow.  Weak  correlations  are  found  in  the  perpendicular  direction. 
These  represent  formation  of  channel-like  features,  with  a  characteristic  size  larger 
than  the  grid  spacing.  The  medium  is  anisotropic,  with  vertical  high  permeability 
channels  spanning  the  box. 


Melt  fraction  decreases  upstream  (conductively  cooled  conditions).  Fluid 
enters  at  z  =  0  with  solute  concentration  c  =  1.  At  the  entering  height  the  fluid  is  in 
equilibrium  with  the  matrix.  The  equilibrium  saturation  decreases  linearly  upstream 
so  that  Ce,(z)  =  1  —  zjL  where  z  is  the  height  and  L  is  the  length  of  the  box.  This 
renders  the  fluid  constantly,  very  slightly,  supersaturated  with  respect  to  its  height.  A 
simple  Ccdculation,  similar  to  that  presented  in  the  calculation  of  the  steady-state  of 
Chapter  III,  can  be  made  to  prove  this.  The  initial  porosity  field  for  this  experiment 
is  not  a  random  distribution,  but  instead  the  “channelized”  porous  media  formed  in  a 
previous  dissolution  experiment.  (The  section  shown  in  Figure  4-5a  is  a  section  from 
the  initial  flow  field  in  this  experiment.) 

Figure  4-5b  shows  a  slice  parallel  to  the  flow  direction  of  the  ZD  flux  field  after 
800  time  steps.  The  initial  channels  of  flow  have  been  completely  wiped  out,  except 
for  remnants  at  z  =  0,  where  the  fluid  enters  at  equilibrium  and  so  no  reaction 
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occurs.  Figure  4-6a  and  b  show  the  correlation  function  in  the  direction  parallel 
and  perpendicular  to  the  flow,  respectively.  Deposition  ha^  completely  destroyed  the 
initial  correlations  (chcinnels  of  flow)  obtained  by  dissolution,  so  as  to  cause  uniform 
eind  diffuse  flow.  In  Figure  4-5b  one  can  also  observe  the  reduction  in  flux  variability 
reminiscent  of  the  results  presented  in  Figure  4-4(bl). 

The  transition  to  Pyroxene  saturation.  To  simulate  the  transition  from  dis¬ 
solution  to  deposition,  fluid  enters  at  z  =  0  with  c  =  0.  The  matrix  is  increasingly 
soluble  until  z  =  Z//2  by  having  Ce,  =  zjL.  At  z  =  L/2,  Ce,  begins  to  decrease  lin¬ 
early  and  Ce,  =  1  —  z/i.  The  z  profile  of  the  equilibrium  saturation  and  the  resulting 
concentration  profile  (averaged  over  x  cind  y,  after  400  time  steps)  are  drawn  in  solid 
line  eind  stars  respectively  in  Figure  4-7a.  Melt  is  undersaturated  (c  <  c^)  for  more 
than  half  of  the  box,  i.e.  concentration  profile  attains  a  maximum  at  a  higher  z  than 
does  the  equilibrium  concentration.  The  vertical  distance  between  the  maximum  in 
concentration  and  the  point  of  inversion  in  equilibrium  concentration  is  determined 
by  the  Da;  The  deviation  from  equilibrium  is  a  result  of  the  competition  between  ad- 
vection  of  downstream  undersaturated  material  and  equilibrating  chemical  reaction. 
As  reaction  rates  increase  (Da  — »  oo),  the  concentration  profile  approaches  that  of 
Ccq,  but  z(mai(c))  will  remain  greater  than  z(max(ce,))  for  any  finite  flow  rates. 

Vertical  sections  of  the  resulting  3D  flow  pattern  are  presented  in  Figure  4-8.  Flux 
channels  form  in  the  dissolving  region  and  extend  somewhat  into  the  precipitation 
region.  Pressure  sections,  on  the  other  hand,  vary  only  along  the  z  direction,  and 
appear  constant  in  i  and  y.  Thus  only  z  profiles  of  pressure  are  shown. 

Figure  4- 7b  shows  average  (over  x,y)  porosity  and  deviation  from  hydrostatic 
pressure  as  a  function  of  z.  Melt  fraction  (porosity,  plotted  in  circles)  increases  to 
half  the  box,  after  which  a  rapid  decrease  in  porosity  occurs  so  as  to  form  a  low 
permeability  cap  over  the  molten  region.  The  peak  in  porosity  corresponds  to  the 
largest  deviation  of  c  from  Ce,,  attained  at  the  peak  of  Ce,. 

The  deviation  of  “mecisured”  pressure  from  hydrostatic  pressure  is  plotted  as  a 
solid  line  in  Figure  4-7b.  Overpressure  attains  a  maximum  at  the  same  height  as 
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does  c{z),  i.e.,  the  point  of  transition  from  net  dissolution  to  net  deposition.  As 
discussed  above,  the  vertical  distance  between  the  two  points,  the  point  where  porosity 
is  maximized  and  the  point  where  overpressure  is  maximized,  is  determined  by  the 
Da,  and  decrea.ses  cis  Da  — >  oo.  In  this  rigid  matrix,  the  difference  in  porosity 
between  the  molten  and  the  crystallized  regions,  and  the  overpressure  associated  with 
the  transition,  continue  to  increase  with  time.  The  consequences  of  such  a  physical 
scenario  will  be  discussed  in  the  conclusions. 


4.6  Summary  and  conclusions 

In  this  study  I  have  presented  a  new  three-dimensional  computer  model  which  ef¬ 
ficiently  models  Darcy-sccile  aspects  of  flow  and  chemical  reaction,  both  deposition 
and  dissolution,  in  a  porous  medium.  Simulations  of  generic  and  specific  process  were 
performed,  and  the  following  conclusions  for  flow  and  reaction  at  high  Da  and  Pe 
numbers,  are  reached: 


4.6.1  General  conclusions. 

Reaction  in  porous  media  at  high  Da  causes  spatial  recirrangement  of  the  porosity 
in  the  matrix  with  relatively  little  effect  on  the  shape  of  the  porosity  histogram. 
Dissolution  results  in  formation  of  alternating  high  and  low  permeability  regions 
parallel  to  the  main  flow  direction,  while  deposition  acts  to  decrease  correlations 
in  porosity  and  diffuses  the  flow  (Figure  4-4).  It  is  possible  that  reactive  flow  is 
accompanied  by  a  change  in  the  relationship  between  porosity  and  permeability  (i.e., 
in  dissolving  systems  global  permeability  will  be  higher  than  predicted  by  a  measured 
average  porosity,  while  in  depositing  systems  global  permeability  will  be  lower  than 
thus  predicted).  However,  in  my  simulations  this  effect  is  not  observed,  and  the  total 
permeability  of  the  system  is  well  predicted  (within  1%)  by  equation  4.9,  using  the 
average  porosity  within  the  whole  box. 
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4.6.2  Relevance  to  melt  migration  in  the  earths  mantle 

The  simulations  in  this  section  followed  closely  the  proposed  different  scenarios  in 
[36].  The  results  of  simulations  agree  with  predictions  in  that  paper. 

Melt  ascending  beneath  midocean  ridges,  on  an  adiabatic  geothermal  gra¬ 
dient.  It  is  currently  believed  that  the  earths  mantle  rapidly  ascends  below  mid¬ 
ocean  ridges  in  near  adiabat  conditions  and  remains  above  the  point  of  pyroxene 
saturation  for  most  of  its  ascent  path,  continuously  increasing  in  liquid  mass  [72,  36]. 
Under  such  conditions,  simulations  show  that  three  dimensional,  elongated  channels 
of  fast  flowing,  undersaturated  melt  form  spontaneously.  The  vertical  channels  span 
the  box  size,  similar  to  the  two  dimensional  linear  predictions  obtained  in  Chapter 
III. 

As  detailed  in  Chapter  III  and  in  [36]  and  [32],  the  spontaneous  formation  of  such 
channels  may  explain  the  observation  of  chemical  disequilibrium  between  midocean 
ridge  basalt  and  surrounding  depleted  peridotite  [28,  27].  Some  kind  of  channels  have 
long  been  assumed  in  order  to  explain  the  geochemical  observations  [79,  22],  but  a 
formation  mechanism  for  channels  in  the  viscously  deformable  parts  of  the  mantle 
was  not  known. 

Melt  ascending  and  cooling  on  a  conductive  geotherm  (intra-plate  volcan- 
ism  and  subduction  arcs).  Beneath  hot  spot  volcanos,  ascending  melt  follows 
cin  adiabatic  path  to  the  base  of  the  tectosphere,  but  will  cool  at  shallower  levels. 
Thermal  models  of  subduction  zones  indicate  that  melt  initially  heats  and  then  cools 
as  it  ascends  through  the  mantle  wedge.  In  both  cases,  PT  diagrams  as  shown  in 
[36],  reveal  that  melts  rising  through  the  upper  mantle  will  first  increase  and  then 
decrease  in  mass. 

I  find  that  when  liquids  increasingly  crystallize,  flow  becomes  diffuse  and  homo¬ 
geneous,  even  if  initially  focused  into  channels.  This  result  is  in  agreement  with  the 
chemical  signature  of  bcisalts  and  mantle  samples  obtained  from  such  regions  ([36] 
cind  references  therein),  which  show  a  closer  approach  to  equilibrium  and  slower  melt 
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extraction  rates,  in  contrcist  with  midocean  ridges. 

Simulations  show  that,  at  least  in  a  rigid  medium,  a  low  porosity  cap  will  form  at 
the  transition  from  dissolution  to  crystallization.  The  sharp  transition  from  a  region 
with  high  melt  fraction  to  a  region  of  low  permeability  causes  rapidly  increasing  over¬ 
pressurization.  Such  regions,  whose  depth  in  the  mantle  it  may  be  possible  to  calculate 
from  thermodynamic  considerations,  will  be  candidates  for  periodic  hydrofracturing, 
due  both  to  the  high  melt  fraction  and  the  overpressurization  [36]. 

A  relevant  question  concerning  the  initiation  of  hydrofractures  in  this  region,  is 
whether  compaction  forces  can  dissipate  the  excess  pressure  faster  than  the  rate  at 
which  it  increases  due  to  a  buildup  of  a  low  permeability  cap.  It  has  been  argued  that 
in  steady-state,  compaction  will  dissipate  any  stress  differences  larger  than  Apgh,  the 
buoyancy  forces  obtained  over  a  compaction  length  h  [81].  This  calculation  is  not  true 
in  a  medium  constantly  changing  due  to  chenoical  reactions.  Compaction,  which  acts 
to  dissipate  the  overpressure,  will  be  competing  against  crystallization,  which  acts 
to  build  up  pressure.  The  feister  of  these  two  processes  will  determine  the  resulting 
overpressure. 

The  above  findings  are  consistent  with  geological  evidence  for  dunite  dissolution 
features  found  in  mantle  structures  in  ophiolites  preserved  from  adiabatic  upwelling 
beneath  an  oceanic  spreading  ridge  [32],  and  fractures  filled  with  products  of  crystal¬ 
lized,  cooled  basalt  found  in  mantle  sections  preserving  structures  formed  within  a 
conductive  geotherm. 

4.6.3  Open  question  and  future  investigations 

The  model  presented  above  assumes  a  constitutive  porosity-permeability  relation 
given  by  equation  4.9.  This  is  the  first  point  for  suggested  improvement,  since  as 
stated  previously,  it  is  possible  to  deposit  at  strategic  small  necks,  with  a  very  small 
effect  on  porosity  but  a  tremendous  effect  on  permeability. 

How  to  improve  4.9  is  not  immediately  clear.  One  way  is  to  use  particle  based 
models,  such  cis  the  Lattice-Gas  model,  where  porosity-permeability  relations  emerge 
physically,  without  a  need  for  analyticcil  predictions.  But  these  microscopic  models 
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axe  limited  in  size.  Macroscopic  modeling  can  be  improved  by  imposing  constitutive 
laws  which  take  into  account  the  changing  situations  in  low  porosity  and  assume  a 
power  n  that  increcises  as  porosity  is  decrecised  in  equation  4.9.  Another  possibility  is 
to  use  only  a  fraction  of  the  porosity  (which  can  be  termed  the  connected  porosity) 
in  the  flow  calculations. 

Reactive  infiltration  instability  in  ZD.  Similar  to  the  2D  reactive  infiltration 
instability  and  as  a  consequence  of  laboratory  experiments  in  ZD  [17],  one  expects 
the  reactive  infiltration  instability  to  form  in  3D,  but  it’s  form  is  yet  to  be  fully 
understood.  In  experiments,  the  shape  of  the  reaction  front  has  been  observed  to 
have  scale-invariant  fingers  under  certain  conditions  [17]. 

Model  simulations  of  the  front  propagating  in  ZD  show  that  the  front  is  unstable 
for  smeill  porosity  perturbations,  but  forming  fingers  lock  on  the  initial  noisy  distribu¬ 
tion  of  the  porosity.  Increased  diffusion  is  expected  to  change  finger  size,  but  I  have 
not  yet  succeeded  in  determining  the  behavior  of  this  instability  in  parameter  space. 
It  is  necessary  to  simulate  such  a  phenomena  on  larger  boxes  in  order  to  observe 
possible  scale  invariant  behavior. 

Melt  flow  in  the  mantle.  The  channeling  instability  in  a  gradient  of  solubility 
may  have  interesting  forms  in  3D,  which  need  to  be  investigated  as  a  function  of  the 
vaxious  parameters.  This  part  of  the  study  is  complicated  in  a  rigid  medium,  because 
of  the  lack  of  a  steady  state,  and  so  future  addition  of  compaction  is  proposed.  Such 
tests  using  a  numerical  code  that  incorporates  compaction,  as  well  as  reaction,  will 
also  investigate  the  trcinsition  from  dissolution  to  precipitation  and  determine  the 
maximum  pressure  obtained  as  a  function  of  competing  rates. 

An  additional  point  to  study  is  temperature  effects.  In  the  earths  mantle,  the 
equilibrium  concentration,  and  so  the  transition  from  dissolution  to  precipitation,  are 
a  function  of  temperature  as  well  as  pressure.  For  example,  fast  flowing  melt  in  a 
high  permeability  channel  may  be  still  hot  enough  to  be  able  to  dissolve  the  solid 
matrix,  while  the  neighboring  regions  may  have  already  cooled  enough  to  precipitate. 
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resulting  in  focusing  of  the  upwelling  melt  in  narrow  chimneys  with  solidified  walls, 
similar  to  results  of  the  experiments  reported  in  [36]  and  [82].  In  the  future,  it  would 
be  interesting,  though  somewhat  difficult  due  to  the  additional  length  and  time  scales 
involved,  to  incorporate  temperature  effects  into  the  existing  model. 
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2d  Linear  stability  theory 

theory  in  solid,  experiments  in  circles 


wavenumber  (in  units  of  l/dx) 


Figure  4-1:  Growth  rate  as  a  function  of  wavenumber  for  the  reactive  infiltration 
instability.  Results  from  model  simulations  in  circles,  analytical  prediction  from  linear 
stability  from  [57]  (equation  4.24)  in  solid  line.  Simulations  were  performed  on  a  box 
with  33^  grid  points.  Parameters  used  are:  ^/  =  2,  =  1, 7  =  0.125,  Vf  =  0.4,  a  =  80 

and  Da  =  40.  The  highest  wavenumbers  correspond  to  features  close  to  the  grid 
spacing.  Such  high  frequency  features  are  not  well  resolved  by  partial  differential 
equation  solvers  such  as  this  model. 
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2.0 


time  (in  units  of  delta  t) 


Figure  4-2:  A  supersaturated  fluid  reacting  and  flowing  in  a  uniform  media.  Equation 
4.25  describes  the  time  evolution  of  the  uniform  concentration  towards  equilibrium. 
Parameter  set  used:  Da  =  10,  A  =  0.63,^  =  l,co  =  0.1,  At  =  0.0125.  Prediction 
from  4.25  in  solid  line.  Simulation  results  in  circles. 
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Figure  4-3:  Average  porosity  as  a  function  of  time  for  (al)  dissolving  porous  media 
and  (bl)  porous  media  undergoing  deposition.  Porosity  in  both  cases  approaches  a 
limiting  value;  in  a  dissolving  matrix  (f)  <f>f  (the  porosity  when  no  more  soluble 
material  is  present),  and  in  a  depositing  porous  media  (f)  0.  Time  is  counted  in 

units  of  initial  time  steps.  (Time  steps  in  model  are  changed  according  to  stability 
requirements,  where  for  larger  fluxes  one  needs  to  take  smaller  steps.)  The  bottom 
half  of  the  figure  shows  2D  sections  (taken  parallel  to  z,  the  main  flow  direction)  of 
the  final  flow  fields,  obtained  in  (a2)  dissolving  and  (b2)  depositing  porous  media, 
lighter  (darker)  shadings  correspond  to  higher  (lower)  than  average  fluxes.  (a2)  shows 
formation  of  elongated  flow  structures  p2irallel  to  flow  direction.  (b2)  shows  more  dif¬ 
fuse  flow  structures,  although  isotropy  exists  because  of  the  driving  pressure  gradient 
cilong  the  z  direction. 
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(a1)  (a2)  (a3) 


Figure  4-4:  Histograms  of  flux  and  porosity,  and  correlation  functions  of  porosity  in 
a  dissolving  and  “clogging”  porous  media.  For  the  histograms  the  mean  value  of  the 
variables  was  subtracted  to  enable  compiirison.  Bin  size  is  constant.  The  correlation 
functions  of  porosity  shown  are  of  lags  along  the  z  direction.  Experiments  performed 
on  a  33^  node  box  with  Da  =  20,  Pe  =  10^  and  vq  =  0.5.  Initial  configuration 
in  solid  line.  Fined  configuration  in  stars,  a.  Dissolving  porous  media:  (al)  The 
variability  in  flux  after  700  time  steps  is  wide  compared  to  the  initial  variability,  as  is 
expected  for  formation  of  focused,  non-uniform  flow  with  regions  of  fast  and  slow  fluid 
velocity.  (a2)  The  change  in  porosity  variability  between  t  =  700  and  t  =  0  is  nearly 
unnoticeble.  This  is  expected  when  permeability  is  changed  due  to  organization  of 
the  porous  medium  and  formation  of  long-range  correlations,  and  not  due  to  changes 
in  porosity  variability.  (a3)  Porosity  correlation  function  in  the  direction  parallel 
to  flow  is  considerably  enhanced,  consistent  with  (al)  and  (a2).  b.  Deposition  in 
porous  media:  (bl)  The  variability  of  the  flux  after  600  time  steps  is  much  smaller 
than  in  the  initial  configuration,  which  means  that  the  flow  becomes  more  uniform 
with  time.  (b2)  The  decrease  in  porosity  variability  is  less  prominent,  as  expected 
when  permeability  is  changed  due  to  recirrangement  of  the  porous  media  to  form  anti¬ 
correlations  (i.e.  preferential  clogging  in  regions  close  to  high  permeability  regions). 
(b3)  Deposition  results  in  reduced  correlations  in  porosity  which  act  to  diffuse  and 
homogenize  the  flow. 
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Figure  4-5:  Vertical  slices  of  the  ZD  flux  field  in  a  reacting  porous  medium  with  a 
gradient  in  solubility.  Box  is  of  33^  grid  points.  Parameters  used:  Da  =  20,  Pe  = 
10^,  Uo  =  0.5.  (a)  Dissolving  porous  medium  with  solubility  increasing  with  z  so  that 
fluid  is  slightly  undersaturated  everywhere.  Initial  porosity  distribution  is  random. 
The  slice  presented  is  taken  after  700  time  steps.  Elongated  channels  are  formed,  with 
a  large  variation  in  flow  velocities,  (b)  Using  the  configuration  obtained  by  dissolution 
as  the  initial  condition,  the  slice  presented  is  taken  from  a  porous  media  after  800 
time  steps  of  deposition.  Flux  channels  have  been  obliterated.  Flow  is  uniform  with  a 
small  Vciriability  around  the  mean.  (Some  trace  of  channels  can  be  observed  at  z  =  0 
where  there  is  no  reaction). 
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Figure  4-6:  Correlation  function  of  the  flux  field  in  a  reacting  porous  medium  with 
a  gradient  in  solubility.  Box  is  of  33^  grid  points.  Parameters  used:  Da  =  20,  Pe  = 
10^,  Uo  =  0-5.  The  initial  condition  for  dissolution  is  a  random  porous  medium.  The 
initieil  condition  for  deposition  is  the  channelized  porous  media  from  the  dissolution 
experiment,  (a)  Correlation  functions  parallel  to  flow  direction  (lag  in  z):  Flux  is 
highly  correlated  in  the  dissolved  matrix,  represented  by  circles,  while  it  is  completely 
uncorrelated  in  the  deposited  porous  media  (in  plusses),  i.e.  the  channels  caused  by 
dissolution  have  been  preferentially  clogged,  (b)  Correlation  functions  perpendicular 
to  flow  direction  (lag  in  x):  In  a  dissolved  porous  media  (in  circles)  flow  is  somewhat 
correlated  in  the  direction  perpendiculcir  to  the  forcing,  i.e.  pipes  have  a  width  larger 
than  the  grid  spacing.  These  correlations  are  seen  to  be  reduced  by  deposition  (in 
plusses). 
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Figure  4-7:  Results  from  simulations  of  the  transition  from  dissolution  to  deposition, 
after  400  time  steps.  Porosity  wets  initially  homogeneous  and  random  with  average 
value  of  1.  Parameters  are  as  in  Figure  4-6.  Plotted  are  various  z  profiles  obtained  by 
averaging  over  x  and  y.  (a)  Equilibrium  concentration  cis  a  function  of  z,  in  solid  line, 
and  resulting  concentration  in  staxs.  the  fluid  is  increasingly  undersaturated  until  a 
certain  point,  after  which  it  is  slightly  supersaturated.  The  vertictd  distance  between 
the  meiximas  in  concentration  and  in  equilibrium  concentration  is  determined  by  the 
Da  number,  (b)  The  porosity  profile,  in  circles,  has  evolved  to  form  a  highly  perme¬ 
able  region  below  a  low  porosity  cap.  Deviation  of  resulting  pressure  from  hydrostatic 
pressure  is  plotted  in  solid  line.  While  porosity  attains  a  maximum  at  the  same  point 
as  the  equilibrium  concentration,  the  overpressure  attains  a  maximum  at  the  tran¬ 
sition  between  dissolution  (undersaturation)  and  deposition  (supersaturation).  The 
deviation  between  these  points  is  thus  also  dependent  on  Da. 
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Figure  4-8:  Two  different  vertical  sections  of  a  2D  flux  field  in  a  transition  simulation. 
Flux  channels  span  the  dissolution  region  and  may  affect  the  crystallizing  region 
as  well.  The  (^  =  1  contour,  which  represents  the  transition  point,  is  drawn  as  a 
horizontal  solid  line. 
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Chapter  5 


Conclusions 


Motivated  by  Vcirious  observations  in  sedimentary  and  igneous  rocks,  this  thesis 
presents  a  study  of  changes  that  occur  in  a  porous  medium  subjected  to  flow  and 
chemical  reaction.  Coupled  flow  and  reaction  operate  during  sedimentary  rock  for¬ 
mation,  during  upwelling  of  melt  in  the  mantle,  at  the  core-mantle  boundary,  and  in 
various  other  technological  and  geological  settings. 

The  coupled  physical  process  influences  the  geometry  of  a  porous  matrix  at  both 
the  microscopic  and  the  macroscopic  levels.  At  the  smallest  scale  flow  may  have  a  role 
of  keeping  the  system  out  of  equilibrium  so  as  to  allow  for  growth  and  existence  of 
non-equilibrium  features.  At  larger  scales  flow  may  act  as  a  “long-range  messenger” , 
causing  Icirge-scale  orgeinization  during  cementation  or  dissolution. 

This  thesis  reaches  two  levels  of  conclusions:  The  first  level  is  that  generic  aspects 
of  flow  and  reaction  change  the  statistical  characteristics  of  a  porous  medium,  with 
dissolution  and  deposition  having  qualitatively  opposite  effects  (Chapter  IV).  It  may 
be  possible  in  the  future  to  use  the  statistics  to  obtain  quantitative  constraints  on 
the  processes  that  different  porous  media  have  undergone.  It  is  also  clear  that  flow 
cind  reaction  effect  permeability,  in  a  way  which  needs  to  be  further  studied.  Future 
studies  require  deeper  understanding  and  quantification  of  the  interaction  between 
the  scales. 

On  the  second  level  this  thesis  suggests  that  coupled  flow  and  reaction  are  sig¬ 
nificant  forces  in  geology:  This  generic  physical  process  may  be  responsible  for  the 
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observed  fractal  structures  on  tbe  pore  sccde  of  sedimentaxy  rocks  (Chapter  II,  also 
in  [1]),  for  formation  of  channels  in  the  mantle  (Chapter  III,  also  in  [3],  and  [36,  32]), 
and  for  controlling  modes  of  melt  extraction  (Chapter  IV  and  [36]).  The  coupled 
process  may  influence  rates  of  lava  flow  and  sea  floor  accretion  cis  well  as  volcanic 
eruptions.  Future  extension  of  these  studies  to  more  realistic  conditions,  both  via 
experiments  and  through  theory,  is  necessary  before  these  new  results  can  be  used  in 
a  qucintitative  fashion. 
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Appendix  A 


Theoretical  prediction  of  the 
roughness  amplitude  (Chapter  II) 


In  this  appendix  we  calculate  a  theoreticiil  value  for  A{D),  the  prefactor  in  equation 
2.17,  for  interfaces  formed  by  Model  I.  Our  calculations  will  use  the  characteristics  of 
the  ’single-step’  model  «ind  the  fact  that  power  spectra  of  self-affine  interfaces  follow 
a  power- law  [87]. 

In  the  calculation  of  A{D)  the  first  step  is  to  define  a  two-dimensional  discrete 
Fourier  trcinsform  H{p,q)  of  the  height  h(m,n),  where  m  and  n  are  the  discrete  x 
and  y  positions  respectively: 

^(P.7)  =  A  E  E  (A.l) 

m=l  n=l 

The  inverse  transform  is 

A(m,  n)  =  E  i;  H(p,  (A.2) 

p=i  g=i 

Since  in  Model  I  the  height  difference  between  adjacent  sites  is  always  constrained 
to  be  of  amplitude  one,  the  sum  of  the  squares  of  the  slopes  in  the  x  direction  is 

E  (A.3) 

m=:l  n=l 

By  substituting  the  inverse  Fourier  transform  for  h  in  the  above  ‘sum-rule’  one  obtains 

2EEI-^(p>7)P(1  -cos(27rp/Iy))  =  1  (A.4) 

p=l  g=l 
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The  same  calculation  in  the  y  direction  gives 


2  E  E  Wp.  «)P(1  -  (2’'3/i))  =  1. 

p=l  g=l 


(A.5) 


and  by  adding  A. 4  cind  A.5,  we  obtain 


S  £  l^(P>  9)I^(2  -  cos  {2'irq/L)  -  cos  i2irp/L))  =  1. 

p=l  9=1 


(A.6) 


On  the  other  hand,  the  width,  cls  defined  by  equation  2.4,  is  also  related  to  the  power 
spectrum  by  the  Parseval- Rayleigh  relation: 


4  E  E  =  EE 

m=l  n=l  p=l  o~l 


p=l  9=1 


Thus  A{D)  can  be  extracted  from  equations  A. 7  and  2.17: 


L-l L-l 


A\D)  =  L^‘>-^)-£'£\H(p,g)\K 


p=l  9=1 


For  self-affine  interfaces  the  power  spectrum  scales  like  a  power-law, 


\H(p,q)\^  =  AXD)  , 


where  A'{D)  is  the  amplitude  of  the  power  spectrum  [87]. 

Equation  A. 8  can  be  rewritten  using  \H{p,q)Y  from  equation  A. 9: 


AW  =  E  E  ((^)'  +  (^) 


2\  D-i 


{A.7) 


(A.8) 


(A.9) 


(A.IO) 


A'{D)  can  be  derived  by  substituting  the  power  spectrum  given  in  A.9  in  A.6  to  give 
a  final  form  for  the  square  of  the  roughness  amplitude 


A^D)  = 


Ej.1  Ej'=i(2  -  cos(2Tg/i)  -  cos(27rp/i))(52  +p=)“-<' 


(A.ll) 


The  sum  can  be  evaluated  numerically,  to  yield  a  theoretical  value  of  A(D),  which  is 
shown  in  Figure  2-15.  Our  calculations  of  A  have  not  fully  converged  and  show  a  very 
weak  dependence  on  L,  with  L  values  taken  up  to  10'*.  The  integral  approximation 
of  the  sum  in  equation  A.ll  shows  however  no  dependence  on  L,  and  therefore  full 
convergence  is  expected  as  L  — >  oo. 
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Appendix  B 

Discussion  of  Parameter  Values 
for  the  Mantle  (Chapter  III) 


Table  3.1  illustrates  the  derivation  of  a  range  of  Damkohler  numbers  which  may¬ 
be  applicable  to  porous  flow  of  melt  in  the  mantle.  The  results  axe  tabulated  in 
two  ways,  in  terms  of  an  “equilibration  length,”  and  in  terms  of  Damkohler  number 
for  a  fixed  system  length  scale  of  100  m.  In  our  study,  the  equilibration  length, 
ieq  =  (f>pfWolR^ff,  is  the  length  over  which  fluid  will  advect  before  equilibrating 
with  its  surroundings. 

Critical  input  parameters  in  determining  equilibration  lengths  for  the  mantle  are 
the  crystal  dissolution  rate  and  the  melt  flow  velocity.  We  present  a  broad  range 
of  possible  values,  because  the  estimation  of  effective  dissolution  rates,  porosity  and 
fluid  velocity  cire  so  uncertain.  Several  studies  have  measured  dissolution  of  mantle 
minerals  in  basaltic  melts  at  upper  mantle  pressures  and  temperatures  [41,  12,  96]. 
All  of  these  indicate  that  measured  dissolution  rates  are  diffusion  controlled,  and 
thus  depend  on  the  diffusivity  of  the  dissolving  species  and  the  width  of  a  chemical 
boundary  layer  around  the  dissolving  crystal.  The  first  two  studies  emphasized  re¬ 
sults  for  relatively  “well-stirred”  melts,  with  narrow  chemical  boundary  layers,  while 
[96]  attempted  to  minimize  convection  and  mixing  in  liquids  surrounding  dissolving 
crystals,  maximizing  the  width  of  the  chemical  boundary  layer.  It  is  difficult  to  know 
which  of  these  apply  to  the  microscopic  geometry  of  melt  flow  in  the  mantle,  below 
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the  continuum  or  Dcircy  scale,  in  which  some  intergranular  pores  could  be  effectively 
stagnant,  while  others  may  carry  rapidly  moving  liquids.  For  this  reason,  we  use 
linear  dissolution  rates  derived  from  these  studies  ranging  from  10~®  to  10"^^  m/s. 
We  convert  volume  to  weight  units  using  an  approximate  density  of  3000  kg  m“^  to 
obtain  the  reaction  rate  constant,  R. 

To  obtain  effective  dissolution  rates  (kg  s“^  ni“^),  we  calculated  effective  surface 
areas  over  which  dissolution  could  occur  in  mantle  peridotites,  per  unit  volume.  This 
calculation  requires  estimation  of  the  grcdn  size,  grain  shape,  proportion  of  solid/liquid 
surface  area  to  toted  surface  area,  and  proportion  of  soluble  phases  in  the  peridotite. 
Pyroxene  is  much  more  soluble  than  olivine  in  ascending  liquids  (e.g.,[33,  36]),  so  for 
the  purposes  of  this  calculation  we  assumed  that  the  proportion  of  soluble  phases  was 
the  proportion  of  pyroxene  in  mantle  peridotite.  We  have  used  the  work  of  [90]  to 
estimate  solid-liquid  surface  areas  for  bascdt-mantle  systems  as  a  function  of  porosity. 
In  calculating  the  solid  surface  ajea  per  unit  volume,  we  have  assumed  cube-shaped 
grains  with  linear  dimensions  from  0.01  to  0.5  cm. 

[36]  Ccdculated  peridotite  solubilities  in  typical  melts  along  likely  adiabatic  P- 
T  gradients  beneath  mido-ocean  ridges,  using  a  thermodynamic  model  for  partially 
molten  silicate  systems.  Thus  this  Ccdculation  incorporates  the  effect  of  the  heat  of 
fusion  in  limiting  solid  solubility.  Results  of  these  calculations  were  used  to  estimate 
an  approximate  value  for  the  solubility  gradient,  given  in  Table  3.1.  This  calculation 
did  not  include  the  possible  effects  of  advective  heat  transport  by  rising  melt  in 
high  permeability  channels.  Potentially,  if  melt  flux  becomes  large  enough,  this  could 
result  in  local  heating  of  channels  to  temperatures  higher  than  the  adiabat  for  partially 
melting  mantle  peridotite.  If  this  occurred,  it  would  increase  the  local  solubility  of 
solid  phases  to  values  higher  than  those  in  Table  3.1,  and  act  to  further  enhance 
growth  of  chemnels. 

Steady  state  melt  flow  velocities  were  calculated  using  the  reasoning  of  [78],  in 
which  the  Darcy  flux,  <j)WQ  =  FVq  (where  F  is  integrated  mass  fraction  of  melting  and 
Vo  is  solid  upwelling  rate)  for  ascending  mantle  beneath  a  spreading  ridge.  Another 
simplification  can  be  made  if  ^  is  constajit  (e.g.,[27,  75]).  In  this  formulation,  if  <f)  is 
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much  smaller  than  F,  flow  velocities  are  much  greater  than  if  the  porosity  is  of  the 
same  order  as  the  melt  fraction.  A  range  of  values  is  used  in  Table  3.1  to  investigate 
the  maximum  and  minimum  equilibration  lengths  likely  in  the  mantle.  Equilibration 
lengths  Ccilculated  in  this  way  remge  from  Angstroms  to  meters. 

The  attempt  to  quantify  effective  reactive  surface  area  brings  about  an  interest¬ 
ing  hypothesis  about  the  finite  size  behavior  of  these  channels:  When  high  poros¬ 
ity  channels  form,  the  solid-liquid  surface  area  per  unit  volume  is  reduced  and  the 
characteristic  equilibration  length  of  the  system  will  increase.  Hence  the  horizontal 
wavelength  of  the  cheinnels  will  grow  (from  equation  3.55).  Increasing  melt  velocity 
with  height  results  in  an  additional  increase  in  Teq-  A  tentative  sketch  of  finite  size 
behavior  is  given  in  Figure  3-9. 
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